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1. Introduction. In the following we consider solutions of the heat conduction 
equation 
oT aT ie. 
at = k ox" ’ k na K, pc, (1.1) 
for zx > 0,7 > O, corresponding to certain heat inputs h(t) for ¢ > 0 over the plane 
zx = 0: initially T vanishes for x > 0. In (1.1) pe is the specific heat per unit volume, K 
the conductivity. 
To this end we start with the “Green’s function” or the “instantaneous heat source” 
solution 
exp [—2°/4kt] 
Gz,)=-) aeko™ °° §?% (1.2) 


0 , t<0. 


The function G satisfies Eq. (1.1) for ¢ > 0 and represents the temperature due to an 
amount of heat discharged at the time ¢ = 0 at x = 0, in a medium of initial temperature 
T = 0, the quantity of heat per unit area of the plane x = 0 being such that 


[ Ge, od@=1, t>0. (1.3) 
The function G is Gaussian in x for each t > 0 and has a deviation varying as ¢'”. 
For « = 0,t > 0, G varies as t”’”. At x = 0, t = 0, G possesses a singularity. 
Assume that in a semi-infinite medium zx > 0, initially at 7 = 0, heat of amount 
h(t) is fed in at x = O for t > 0. The temperature is given for t > 0 by the following 
definite integral: 


Ke, 6 = = [ n(eye(e, t — ¢) av 
pec 70 (1.4) 





_ 1 f' pn exp [—2'/4kt — ¢)] 2, 

aa 4 | h(t’) [k(t ae t’)|'”? dt ° 
The factor pe in (1.4) is due to the specific heat of the material per unit volume: the 
factor 2 in the first integral is due to the fact that in (1.2) the heat flows to both sides of 
x = 0, while h(¢) is defined as the heat flowing only to the side z > 0. 


*Received February 23, 1959. 
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For x = 0, Eq. (1.4) yields 


ae ay @ h(t’) dt’ 
7 0, a SSDS | . 2:7;  «4ax«rkxsise” (1.5 
©, i) ss ixk(t — t’))'? (1.5) 
In particular, let 
(*/n! = f* _- 
h(t) Fi n! t"/T(n + 1), t>.0, n>-—l. (1.6) 
0, t< 0. 
Upon introducing the variable of integration 
t’ = tu, (1.7) 
Eq. (1.6) may be reduced to the beta-integral, yielding 
"T(1/2) "Kael 
l/2 Raienriseouieialls (1.8) 


T0,) = — SS ES - 
ie "Tin + 3/2) ~ (peKr)™ (1/2) - ‘aI 2)(n + 1/2 


Equations (1.6), (1.8) are valid even for fractional n > —1/2, provided n! is in- 
terpreted as T (n + 1). 

The explicit expression (1.8) for T(0, é) can be applied to general h(t) by approxi- 
mating to the latter by means of a polynomial in ¢ 


h(t) = ho + hit + hot?/2! + --- +h,t'/n! (1.9) 


and carrying out the corresponding superposition of the solutions (1.8) 


pi/2 ho hyt No ol 
T(0, t) = atom as 1/2) * 41/2)(3/2) * 0/2)3/2 








1.10) 





-(n + 1/2) 


Turning to the integration of (1.4) for general x, it is shown in Sec. 2 that for the 
heat input (1.6) form = 0, 1, 2, --- the resulting temperature is given by 





ne = “pies : x 
Tl =T,(z, t) = (oc) 7? f,(u), u = Xhkt) , (1.11) 
TABLE I 
n Fr. G, 
—] 1/2 0 
0 1 —2u 
] 2(1 + u?)/3 —2u — 4u3/3 


2 (4 + Qu? + 2u*)/15 —u + 4u*/3 + 405/15 
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where f, , as indicated, depends only on u, and is given by 
f,.(u) = 2n7'?P,(u) exp (—u’) + Q,(u) erfe(u) , (1.12) 


where P,(u), Q,(u)/u are certain polynomials in u’ of degree n, and “erfc” denotes the 
“eomplementary error function.” For n = 0, 1, 2, P, , Q, are given in Table I. The row 
n = — 1 in Table I is explained in Sec. 2, where recurrence equations for P, , Q, are 
also given, as well as expansions for 7’, in powers of x. 

Solutions for 7, for non-integer n are discussed in Sec. 3, where operational expressions 
for T,, are also given. It is shown that these solutions of (1.1) can be extended to x < 0 
and correspond to proper initial temperatures which vanish for x > 0. 

2. Solutions for polynomial power inputs. We consider heat inputs at x = 0, of the 
form (1.6) for integer n 


h(t) = h(t) = 1t"/n}, (2.1) 


where 1 = H(t) is the Heaviside unit function defined by 


0, 
\(t) = H(t) = ‘ me es (2.2) 
0 for ¢ <0. 
[t will be noted that h, satisfy the relations 
dh,(t) _ { 92 
-— h,-,(0). (2.3) 
Therefore the corresponding temperatures 7,,(7, ¢) will satisfy similar relations 
OTA, 7 (2,0, Tir, =0 for t<0. (2.4) 


0 t . 


The sequence h,(¢), 7',(¢) may be extended by means of (2.3), (2.4), but not directly 
by means of (2.1), ton = —1, yielding 











dhjt)  dH(t) “ 
h_,( = 7 = 6%), (2.5) 
Ht) 
/+ 
| 
| 
I 
| 
| 
| 
| 
' 
: 
ai 
= 4] € - 
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where 4(¢) denotes the “unit impulse function”, or the “Dirac function”. Indeed, if 
H(t) be approximated by means of an analytic curve as in Fig. 2-1, then its slope will 
take on the appearance shown in Fig. 2-2, showing a hump of unit area near t = 0. 
In the limit, as « — 0, there results an instantaneous heat input for which the temper- 
ature is given, except for a factor 2/pc, by Eq. (1.2), namely, 








Sit) = 
dt 
-€ O a t 
Fic. 2-2 
-1/2 , 
———73 exp (—u)’ “u= _—— P t>0 
9 r\172 © ’ WP , 
T.,(z, 2) = —G(z, = al ual ane (2.6) 
si lo, t<0. 


It will be noted that Eq. (2.6) agrees with Eqs. (1.11), (1.12), provided P, , Q, are 
chosen as in Table I forn = — 1. 
For n = 0, when (2.1) yields h(t) = 1 for ¢ > 0, the temperature 7, may be calculated 


from (2.4), (2.6) or from (1.4). Carrying out the integration by parts, one obtains 
, t \'”) 2 exp (—u’) | 
fort) = (LY [2222 — mates] a 
T(x, t) al | 73 u erfe u 7) 
= 9 


The corresponding integrations (1.4) or (2.4) have been carried out for 7, for n = 1, 2. 
The results suggest for general integer n, the form (1.11), (1.12). Indeed, substitution 
of (1.11) in (2.4) verifies the assumption (1.11) provided the recurrence equations 

[(n + 1/2)f.(u) — (u/2)fn’(u)] = fr-i() (2.8) 


(2n+2) 





are satisfied. Multiplying both sides by 2/u , there results 


(2n thee — f(W) _ 2frs(u) (2.9) 


y 2n+1 2n+2 


’ u U 
where the left side is the derivative of —f,(u)/u’"*’. Hence, 


fu(u) = —2ue"*? | foal ee + Cu, (2.10) 
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where C is a constant. In view of the condition T, — 0 for z — © or t— 0, the choices 
= + o,C = 0, are proper. One obtains 


fu) = Qui” a i fr pale (2.11) 
For n = —1 Eq. (2.6) yields 
f_s(u) = exp (—w’)/x’”. (2.12) 
Hence, Eq. (2.11) now leads to 
ro) 2 
fo(u) = 2", [ exp ( _F du. (2.13) 


Integration by parts again leads to (2.7). 
Equation (2.7) is of the form (1.11), (1.12) with 


P. = i. Qo = —2Qu. (2.14) 
Applying (2.11) for n = 1 yields f, of the form (1.12) with 
P,=204+0), Q = —u— Se’. (2.15) 
3 3 
A similar calculation for n = 2 shows that T, is given by (1.11), (1.12) with 
= 4 3 : 2 4 -_ — 
P, ig +e + yt, 2. = -[u+4 ue + a (2.16) 
Equations (2.14)—(2.16) are summarized in Table I. As pointed out, forn = — 1, 


Eq. (2.6) still agrees with (1.11), (1.12), (2.12), provided we choose P_, , Q_,; as in 


Table I. 
For general integer n, there results upon substituting (1.11), (1.12) in (2.8) the 


following recurrence equations for P,(u), Q,(u) 
(n + 1/2)Q,(u) — uQ,’(u)/2 = Q,-,(u), (2.17) 
(n + 1/2)P,(u) — uPi(u)/2 + wP,(u) + uQ,(u)/2 = P,-,(1). (2.18) 


Equation (2.17) determines Q, except for the term u’"**. Equation (2.18) then de- 


termines this term and P, . 
The relation (2.4) may be applied to express T', as power series in x, by starting with 
the expansion for 7'_, obtained from Eq. (2.6) 


¢7/2 | 4 ut 
= ie l—uwu +¥-..] 
(rpcK) 2! (2.19) 


1 J z x % 3/2 at °/? 
= —4773 le — — _— | 
(xpcK) 4k * 214k) 
and integrating (n + 1) times termwise with respect to ¢. A single integration yields 


2ny-—n+1/2 


bz ae ; bad oo. r t 
To _ r, “ z. k"2?"(n ame 1/2) + go(x), (2.20) 
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where gp is the constant of integration which may depend on x. This may be determined 
by noting that g(x) must satisfy Eq. (1.1), since 7, and the series in (2.20) satisfy it. 
Hence g(x) reduces to a first degree polynomial in x whose coefficients may be de- 


termined from the heat input condition at x = 0 
er 
— : - = h(t), (2.21) 
GL \eud 
and from 
T(0.t 2 > 99 
740, 2) = -> (2.22) 
\ ° (pcKr)' . 
which follows from (1.8) for n = 0. There results 
“he — (—1)" x" xr 
T,.=-—S>T —~——_____ — —.. (2.23) 
> (pcK r)'”” Do mn — 1/2) K 


Further ¢-integrations of (2.23) and similar determination of the constants of inte- 


gration yield 


¢°/2 o (— 1)"u?" 1 ( z) 
= - ae ; aes a ——\|\2 2.24) 
qr (pcKm)'”” n(n — 1 /2)(n — 32) K n+ 3!k (3-24 
; "atid _ (—1)"*'u*" 1 (= xe x ) . 
T. = S773 _ ar - = rom ie (2.25) 
I (pcKmr)’" “4 n'(n — 1/2)(n — 3/2)(n — 5/2) kK \2! + 3k + otk” : 


Similarly, there results for 7,, for any integer n, the termwise integrated series along 


with the polynomial 


1 | xt” 4 a’t"' ’ pe (2.96) 
= a 2 re no he 2.20 
K Ln! (n + 1)tk (Qn + 1)tk 
A further relation of interest between 7, 
9°T (x, t) Eo ™ 
: ae = —7 (2, i (2.27) 
Ox” & 


follows from (2.4) and the fact that 7, is a solution\of (1.1). 
The above expansions, while convergent for all u, converge slowly for large u, hence 
small ¢. For large wu it is preferable to use asymptotic series of the form 


a 12.8) 3 | A A, 
35 = a exp(—v)| soastant+::: dd, (2.28) 
t’“(pcK) T u u 
where A, , Az - are constants. Indeed, for n = —1 such a series follows from (2.6) and 
P ] a l 1.3 
eric (9) = <5 ep (—9) 1 - — a3 tT are ol (2.29) 
T v 2v Zu 
For n = 1, 2, --+ one assumes (2.28) and applies (2.4), (2.29) to determine the co- 
efficients A, , A2, -:: for successive n. 


Direct substitution of (1.11) in (1.1) shows that f,(w) is a solution of the differential 
equation 
f'(u) + Qufi(u) — (An — 2)f,(u) = 0 (2.30) 


which vanishes for u = + o~. 
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3. Half integer and fractional power inputs. Operational expressions. By differentiat- 

ing 7’, with respect to x, one obtains solutions of (1.1) corresponding to the heat input 

(1.7) for values of n differing from integers by one-half. Indeed, consider the function* 

OT (x, t) 

Ms —-—S 3.1 

os” tea 
where 7, with integer ” are as in Secs. 1, 2. The heat input of 7/ at x = 0 is given by 


__xti| _ x PTA, b | 
a) a es ssi 


Since 7', satisfies (1.1), a°7,,/dx° may be replaced by (1/k) (0T,,/dt), and hence, upon 
recalling (1.8), 


_Kat.) __va/aer" . 
se a k dt |,20  (ak)'?T(n + 1/2) (34) 








This proves the above statement regarding T% . 
Recalling the form (1.11) for T,, , one obtains from (3.1) 


Ti = —(t'/2K)fi(u), u = 2/2kt)'”, (3.4) 


n 


and this can also be put in a form similar to (1.12). 
Of special interest is the case n = 0 for which Eqs. (2.13), (3.4) yield 
T! = _ folw) 2 1 = = 
- 2K K 2(kt)'”” 
For x = 0,¢ > 0, this reduces to 1/K. Hence, the function KT{ corresponds to a sudden 
temperature rise at x = 0, equal to 1. As shown in Fig. 3-1, at various instants the 
abscissas are changed in a fixed ratio. The heat input at z = 0 varies as ¢'”’. 
It is of interest to note that K7% can be obtained by dispensing with heat sources, 
= — o and starting with the initial temperature 
( 
= 2H(-1) = { 
(2 for « <0, 


(3.5) 


but extending the medium to x 


T(x, 0 Pe RES (3.6) 


for the broken-line extensions). 


(see Fig. 3-1] f 
27), (3.1) follows that the sequence of functions 7’ , T,, can be similarly 


From (2. 
extended to x < 0. In particular, the functions 


uo = AT<Z/2, u, = KT./2, uw. = EKT; /2, (3.7) 


u 


aa kKT;, -« a = kKT}, 2, 


3 


form a sequence of solutions of (1.1) satisfying the recurrence equations 
OU», OUm 
= —Un-1 ~— 


ce iy Pm Rates (3.8) 


and such that for ¢ = 0 
u,(x, 0) = {* ir 2 >, (3.9) 
(—x)"/m! for «<0. 


*The — sign is chosen to render 7°/ positive. 
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This may be verified directly by considering (1.11), (1.12) for negative zx, letting t — 0, 
u—— o, and noting that efre (— ~) = 2. 
As indicated in (3.7), and noting (3.1), 
Uea+i = Vit. 2; 
un, = —k" K(AT,,/a2)/2. wa 


We may solve for u,, in terms of its initial values (3.9) by means of G from (1.2) 


u,(z, t) = [ Un(s, 0)G(a — 8, t) ds 
. (3.11) 


0 


== / s” exp [—(x — s)*/4kt] ds/2(rkt)'m!. 


Upon putting s = xz — v(4kt)'”’, there results 








Q"-"(kt)*”? 
Un(z, t) = ——=a7z | Im(u), (3.12) 
where 
u = 2/(4ki)'”, Im(u) = [ (v — u)” exp (—v’) dv. (3.13) 


By expanding (u — v)” by the binomial theorem and integrating v* exp (— v”) dv by 
parts one can again express wu, in terms of exp (— u’), erfc u, and polynomials in wu. 
With the possible exception of the application of the binomial theorem and Eq. 
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(3.10), Eqs. (3.8)-(3.13) apply equally well to all realm > —1, provided m! be interpreted 


as I (m + 1). 
For x = 0 there results (for both integer and non-integer m) 


u,(O, t) = 2"-(kt)"?x/2g,(0)/m!, 
(3.14) 
Mit t) = 2”-*(kt) rng /292 (0) /ml. 


It is of interest to obtain an operational representation of the above solutions. To 
this end we replace the operator 0/d¢ by the symbol p in (1.1): 


aru 


Solving as if p were a constant, one obtains 
u = B exp [2(p/k)'”*] + A exp [—2(p/k)'”], (3.16) 


and dropping the first exponential for x > 0 (presumably because otherwise u = © for 
r= +), 


u = A exp [—(p/k)'”z]. (3.17) 

This yields for z = 0 
u(0, t) = o(t) = A, (3.18) 
_-K ou _ =O = Ka(®) = (Kpc)'’(p)'7A. (3.19) 

Hence 

h(t) = (Kpe)'*p*"o(t), (3.20) 
v(t) = (Kpe)'*p'*h(t). (3.21) 

Interpreting p” as 
p h(t) = es [ (t — 8)""h(s) ds, (3.22) 


for both integer and non-integer n > 0, there results 


ail — 1 * _ h(s) , 
v(t) = T(0, t) Pm (K pe)” nie L (t pas s)'” ds (3.23) 





which agrees with Eq. (1.5). 
Putting (3.20) in the form 


h(t) = (Kpc)'*plp"'*0()], (3.24) 


there results 


1 > 12 a ‘wos 
NO = Faq Ko)? 5 [ Gam ne ds. (3.25) 
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Equation (3.23) is an Abel integral equation for h(t), and Eq. (3.25) yields its solution 
(see for instance, [1]). 
For heat input given by 


h(t) = ¢/Tin+ 1) =p'"l, (3.26) 


where n is either integer or fractional, Eqs. (3.17), (3.19) yield 


: ] 1/2 ] "tinea 
(0, 2) =o) = : ‘ = — 5 Seo na ‘ ijo.ee 
T.,.\ v(t) (Kec)? P 1 (K pe)? Tin + 1/2) (3.27) 
while Eq. (3.17) yields 
as ! _ exp [—(p k)' 2] — 
i. = (Koc) Tin + 172) pri E. (3.28) 


By interpreting these operational expressions as Brownwich integrals in accordance 
with 


. 


l * i ae 
f(p)l = 5 I 5 MP) a, (3.29) 


where L is a proper path of integration in the complex p-plane, one obtains an alternative 
(contour) integral representation for the solutions 7,, and hence also for u,, . 
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RESPONSE OF SHALLOW VISCOELASTIC SPHERICAL SHELLS TO 
TIME-DEPENDENT AXISYMMETRIC LOADS* 


By 
P. M. NAGHDI, University of California, Berkeley, California 
AND 


W. C. ORTHWEIN, International Business Machines Corp., Oswego, N. Y. 


1. Introduction. The quasi-static treatment of problems in the linear theory of 
viscoelasticity has received increasing attention in recent years. The method of solution 
employed in such problems rests on the use of the Laplace transform (to eliminate the 
dependence on time), and the correspondence principle—between the field equations 
and boundary conditions in the linear theories of homogeneous and isotropic elasticity 
and viscoelasticity—which, in the absence of thermal effects, has been established for 
incompressible media by Alfrey [1], and in general form by Lee [2]. The extension of 
Lee’s analogy to problems involving time-dependent temperature fields has been very 
recently given by Sternberg [3]. Also, considerable attention has been given to oscillation 
and wave propagation problems of viscoelasticity in which the inertia terms have been 
included, e.g., [4, 5, 6], and additional references on the subject may be found in a recent 
survey by Lee [7]. 

Closely related to the scope of the present investigation is the recent work on vibra- 
tions of thin shallow elastic shells by E. Reissner [8], who, by utilizing the linear differ- 
ential equations due to Marguerre [9], has shown that for transverse vibrations of shallow 
shells the longitudinal inertia terms (with negligible error) may be omitted; and hence, 
the formulation of the elastokinetic problems of shallow shells, as in the case of elasto- 
staties, may be reduced to the determination of axial (or transverse) displacement and 
an Airy stress function. Subsequently, E. Reissner [10] dealt with transverse vibrations 
of axisymmetric shallow elastic spherical shells, and in particular, obtained the solution 
for an unlimited shell due to an oscillating point load (varying harmonically in time) 
at the apex. 

The present paper is concerned with the response of shallow viscoelastic spherical 
shells to arbitrary time-dependent axisymmetric loads; the medium is assumed homo- 
geneous and isotropic. Although emphasis is placed on unlimited shallow spherical shells, 
shallow spherical shell segments are also considered and discussed in Sec. 7. The solutions, 
employing the differential equations governing the transverse motion of thin shallow 
elastic shells, are obtained with the joint use of the Laplace and the Hankel transforms, 
which, by interchanging the order of the inversions, avoids an otherwise intricate task 
of contour integration in the complex Laplace transform-plane. Explicit results in integral 
form are deduced for viscoelastic shells under instantaneous pulse loading (including 
those uniformly distributed about and concentrated at the apex), and are particularized 
to the cases of Maxwell and Kelvin solids. The solution for a shallow elastie shell and 


*Received November 28, 1958; revised manuscript received May 12, 1959. The results presented 
in this paper were obtained in the course of research sponsored by the Air Force Office of Scientific 
Research under Contract AF 18(603)-47, when both authors were at the University of Michigan. 
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for the case of a flat plate are also given as by-products of the general solution and com- 
parison is made with known results [10]. It may be further noted that the transform 
technique employed here appears to be useful also in connection with other axisymmetric 
problems of stress wave propagation in viscoelastic solids. 

2. Preliminary background. With reference to rectangular Cartesian coordinates 
x; , the stress-strain law for an isotropic and homogeneous viscoelastic medium may be 


written as’ 
lay —_— vi 
P,(0)8;; rem E 2( Ae; ; (2.1) 
P;(O)o3; = P()¢;; ’ 


where o;; and ¢,; are the components of the stress and the strain tensors s,; and e,; 
designate the deviatoric components of stress and strain, the operators P,,(@) involving 
the constant coefficients C“”’ (m = 1, 2, 3, 4) are defined by 


Nm 
Cn 6"; [Cn 0) 


P,,(0) = 
wis (2.2) 
ee 
ilies at” ’ 


and #é denotes time. 
For future reference, we also recall that the Laplace transform with respect to t of 


a (suitably restricted) function U(z, t) is given by’ 
U(x, 8) = L{U(z, #); 8} =| e“"'U(a, t) dt, (2.3) 
0 


where s is the transform parameter, and that the Hankel transform of order zero of the 
function U’(zx, s) is defined by*® 


ure,» = f ” eJo(gx)U"(a, 8) de (2.4) 


provided that (i) the integral {4 U’(z, s)dz is absolutely convergent, and (ii) the func- 
tion U’ is of bounded variation over the region of interest. Furthermore, in connection 
with the Hankel transform of dU’(z, s)/dx, we need the property that (iii) zU’(z, s) 
vanishes both at z = 0, and as zr > ; this ensures the existence of the inverse transform 
of U*. Here it may be noted that the Hankel transform defined by (2.4) formally differs 
from that defined in [17]; however, with the transformation U’ (z, s) = z~'” V (V being 
the function corresponding to U’ in [17]), after multiplying both sides of (2.4) by ¢'” 
and setting V* = ¢'” U*, the two definitions are brought together. 

Since in the theory of shells (and plates) the stress resultants and the stress couples 
are defined in terms of components of stress, rather than components of stress deviation, 
it is expedient for our present purposes to obtain an alternative form of (2.1). To this 


1The Latin indices have, unless otherwise stated, the range of i, j, = 1, 2, 3, and the repeated in- 
dices imply the summation convention. 
_ %See, for example, Churchill [11]; the argument z in U refers to the space variable. 
*See, for example, Sneddon [12]. 


1960] SHALLOW VISCOELASTIC SPHERICAL SHELLS 109 


end, observing that the operators P;’(6@) are in general noncommutative [13], we take 
the Laplace transform of (2.1) and express o/; in terms of ¢/; by* 


oi; = P.(s)P;"(s)et; — $[P2(8)Pi'(8) — Pals)Ps'(s) Jebe 5:3 , (2.5) 
oi; = P,(s)Ps'(8)ei: , 


where s,; is the Kronecker delta. 

It follows from the correspondence principle that the field equations and the boundary 
conditions governing the original viscoelastic problem are reducible to the field equations 
and boundary conditions of an associated problem in the linear theory of elasticity, 
with Young’s modulus £ and Poisson’s ratio v of the elastic solid replaced by 


et = [P,(s)P,(s) + aP(OPOT"Y SPA(s)P (6) } (2.6) 
v(s) [P1(s)P. as) — P 2(8)P 3(8)] 


It is easily verified that the correspondence principle and the results (2.6) are also valid 
for any of the various (consistent) theories of thin shells. 

We also note that in (2.6) the linear, homogeneous and isotropic elastic medium 
may be identified by allowing E(s) — E and v(s) — »v (corresponding to P,(s) = P,(s) = 
i, P;(s) = 2u, P,(s) = 3K, » and K being the shear and the bulk moduli of the elastic 
solid, respectively). The operators P,,(s), when associated with the Maxwell solid, are 








— of, = 
reste. ee ae (2.72) 
P; = 8, P, = 3Ks 
and (2.6) becomes 
8 

E(s) = =1 E, 
s+ #1+»)r (2.7b) 
8 + 1 = v et 





Ae) = s+ 21 +y)7" om 
where r = »/u is the relaxation time, 7 being the viscosity. Similarly, for the Kelvin 
solid, 

















P, = 1, P, = 2yu(1 + 7s) (2.8a) 
P; = 78, P, = 3Krs 
and 
E(s) = 1 + 78 2 
1+ i-2 TS 
3 
1 — (2.8b) 
1— TS 
3v 
vs) = es we 
1+ Ts 


Ww where i in (2.8b) r denotes the retardation time. 


‘It should be noted that unlike P;'(@) the operators P;1(s) are commutative. 
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3. Differential equations for transverse vibrations of shallow elastic spherical shells. 
Let H denote the rise of the shell segment, L the characteristic length (which for spherical 
shells may be conveniently taken as the radius of curvature 2), h the shell thickness, 
p the density, and T the representative wave length. Then, for transverse vibrations 
of thin shallow elastic shells, with the stipulation that (H,/L)° < 1, it has been shown 
by E. Reissner [8] that the effect of longitudinal inertia may be omitted (with negligible 
error) from the differential equations governing the motion of the shell as long as (I'/L) 


r 


and when T is characterized by the following classification 


is of order unity’, i.e., 


| 
(a) If (2 = 0(1) orsmaller, then 


' l ; 
I" = a x (hy), (3.1b) 


, HY? (y? . 
, = 1 (4)(%), 9-1¢) 


where y ~ = pw /E, w being the circular frequency. 

Thus, with reference to cylindrical polar coordinates (r designating the polar radius) 
and with omission of the longitudinal inertia term, the differential equations for the 
axisymmetric transverse vibrations of shallow elastic spherical shells are characterized 


by 110] 


DV?Z72w + 7° 2 = gh ——-— Dr, t 
R p ot f 
Ps hk . 
V'V'F — = Vw = 0, —_ 
ii 
and the various stress resultants and stress couples are given by 
— » op lop 3.3 
\ Si OF ar. \ ‘ / F _ oF or >.) 
/ r 
| . i . 1. ; ow Q « 
M, = —D| 0w/or + Ow / or ; M, = -—D ow/or i , q J wate 
; r or 
0 ' 
Q= —-D—(V'w) 3.3¢) 
or 


where w is the axial displacement, F is the Airy stress function, p is the axial component 


of the surface load, 


5Actually in [8], Reissner has further concluded that the neglect of longitudinal inertia terms is 


justified when (IT'/L) is of order of unity or smaller, but not when (T/L) > 1. 
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a 1a 
ate? >t) 


Eh’ ha 
; and VYV()= > ie 


121 — rv)’ 

It is also relevant to recall here that the steady state solution (for axisymmetric 
transverse vibrations of shallow elastic spherical shells) given by E. Reissner [10], where 
w and F are assumed to have the form 


D = 


w = W(r) exp (tw), F = f(r) exp (iwt), fi = (-1)'”], 


involves Bessel functions (Jo , Yo , Jo , Ko) of argument Ar, where 


= ets (“) | — (Elo) 
ies (Rh)? E ¥ , 7 = 7 (3.4) 


Indeed, since by (3.4) the case of y/R = 1 (corresponding to \ = 0) is not admissible, 
for the axisymmetric vibrations treated in [10] two sets of solutions of (3.2) associated 
with the two frequency ranges R/y 2 1 exist. 

The integration of the second of (3.2), together with the condition of vanishing 
circumferential displacement (which, as in the elastostatic solution of shallow spherical 
shells [14], demands the vanishing of the coefficient of the logarithmic term) leads to 


_ RE ' 7 

Vir = R w - G,(t). (3.5) 

If attention is confined to unlimited shallow shells, then since w(o, ¢) = O and 
VY F(@«, t) = 0 (the latter condition, by (3.3a), is due to independent vanishing of 


N, and N, atr = @), it follows that in (3.5) G,(t) = 0 and for unlimited shallow shells 
the system of differential equations (3.2) reduces to 


DV’?T*w + : V’°F = —ph a + p(r, t) 
) (3.6) 
yp ME 
V'F = Rw 


Furthermore, for unlimited spherical shells, the remaining boundary conditions 
associated with (3.6) are 
Ow 2 oy 
~ (o,t) = Vw(~, t) = 0, (3.7) 
c 


the regularity requirements for oscillating distributed load are specified by 


w(0, t), <0, ), N,0,8, NO, 2d); finite (3.8a) 
M,(0, t), M,(0, t); finite, (3.8b) 


and those appropriate for an oscillating point load are given by (3.8a). 

4. Unlimited shallow viscoelastic spherical shells. The differential equations and 
the boundary conditions (as well as the regularity requirements) appropriate for an 
unlimited shallow viscoelastic spherical shell (O < r < @), subjected to an arbitrary 
time-dependent axisymmetric load, following the application of the Laplace transform, 
and with an appeal to the correspondence principle (Sec. 2), are obtained from (3.6), 











112 P.M. NAGHDI AND W. C. ORTHWEIN (Vol. XVIII, No.2 


(3.7), and (3.8) with the moduli £ and » (of the elastic solid) replaced by E(s) and v(s), 
respectively. In particular, with zero initial conditions, 1.e., 


w(r, 0) = - (r,0) = F(r, 0) = 0, (4.1) 
the differential equations of motion in the Laplace transform-plane read 


D(s)V7?7?w’ + é V°F’ + phs’w’ = p’ 
nB\s) (4.2) 
R ; 

As the solution of the system of differential equations (4.2) involves Kelvin functions 
whose arguments are polynomials in s, thus prohibiting simple inversions, the determi- 
nation of w and F will in general require cumbersome contour integration in the com- 
plex Laplace transform-plane. To overcome this difficulty, we consider the application 
of Hankel transform of order zero to (4.2), and require that w and its derivatives up to 
the fourth, and F’ and its derivatives up to the first vanish to a suitable order at infinity 
such that all integrals employed in the following analysis exist. 

Recalling the formula for the derivative of Hankel transform of order zero [12, 
p. 62], i.e., 


V'F’ = 


[ (V*w’)rJ (rt) dr = —Pw*(é, s) (4.3a) 
and by iteration 

[ (V?V?w')rJ o(ré) dr = E*w*(E, 8), (4.3b) 

70 


then, with the aid of (4.3) and application of the Hankel transform of zero order to 
(4.2), we reach 





D(s)é*w* + ? (V°F)* + phs'w* = p* 
(V'°F)* = va w*, (4.4) 
Elimination of (V’F)* from (4.4) results in 
w* = D“(s)[é* + da(s)]'p*€é, 8), (4.5) 
where 
6 _ 12/1 aol (ze) _ p+ 4 1211-9] 2 
\o = (Rh)? 1+R E@/ | = "+ PEs) ps, (4.6a) 
J 2 
— (4.6b) 


~ 12[1-— »%(s)] 


It is now clear that since the right-hand side of (4.5) involves only polynomials in 
¢ and S, the inversion of w* in s is easily carried out, and its inversion in the Hankel 
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transform parameter £ (being a real variable) is also possible. Thus, by first taking the 
inverse Laplace transform of (4.5), followed by the inverse Hankel transform (which 
is a self-reciprocating transform), we obtain 


w(r, t) = [ EJ o(rt) dé [ L™ {p*(é, 8); S}L™'{y; t — £} dg, (4.7a) 


where 
¥ = D'“(s)[e* + Ads)’. (4.7b) 


It remains to determine the function F. However, since F does not necessarily con- 
form to the requirements for the validity of its Hankel transform, while N, and N, as 
given by (3.3a) are independently finite at r = 0 and vanish at r = ©, we proceed 
instead to establish V’F and 1/r dF /dr. To this end, we turn to the second of (4.4) and, 
with the aid of (4.5), write 


(vin = Ee p{ Die + phe + BON (4.8) 





Again, taking the inverse Laplace transform of (4.8), followed by the inverse Hankel 
transform, we obtain 


- h , . _ -1 + 
VFO.) =e] Bole) ae [ Lip, 8); IL“) VE, 9); ¢— shat, 49a) 


~ ie 
and by integration, since N,(0, t) = N,(0, t) = 3 V’°F(0, t), 
. 1 oF 8a 277 
=s=—=— = = Ot 
i, oa @ 3 ; rV°F (zx, t) dx (4.9b) 


which completes the desired solution. 

Before closing this section, we record some special types of loading, as well as their 
Hankel transforms, which will be of interest presently. For a pulse instantaneously 
applied and removed at ¢ = ¢,, 

P(r, t) = g(r) &(t — &), (4.10) 
where the Dirac delta function is defined by 
ht—t) = 0, (t ¥ tb) 


If, on the other hand, the pulse applied at ¢ = é continues to act indefinitely, then it 
is only necessary to replace 6(¢ — t)) in (4.10) with the Heaviside step function H(t — to), 
defined by 


H(t — t) = 3% (t= b). (4.12)° 
0 (t < to) 
’The use of the letter H with argument (¢ — éo) for the Heaviside step function should not be confused 


with the notation for the shell rise in Sec. 2. 
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Among the various forms which the function g(r) may assume, we specifically cite those 
associated with (a) uniformly distributed pulse over 0 < r < a (a/R < 1), and (b) 


point load pulse applied at r = 0, each given respectively by 
qr) = —pH(a —r), $.13a) 
7 
qr) = —— Or). (4.13b) 
Tr 


Also, for future reference, the Hankel transforms of q specified by (4.13a) and (4.13b) 
are given respectively by [12, p. 88] 


: ofl 
g*(é) = —Fe J ,(aé), (4.14a) 
and by [15, p, 67] 
r 
we = 1.14b) 


5. The elastic solution as a by-product. Reduction to known results. By allowing 
E(s) — E and »(s) — »v in (4.5) and (4.8), and by taking their inverse transformations, 
we deduce the complete solution for the unlimited shallow elastic spherical shell under 
arbitrary time-dependent axisymmetric load. In particular, when the load is specified 
by (4.10) and (4.13a), then with the aid of (4.14a) and tables of integral transforms 
[16, Table 5.2] and [9, p. 323], there follows 








w(r, t) = 7 | Solré) J. (ade aé, (5.1) 
ph Jo 
V’F(r, t) | | pelle | 
>= — me [ J (rt) J (at) dé, (5.2) 
OF , | l [ pak on 
petal - bh ze 1 . 
ar 1 9 | Ce ae” at 
where 
: 3 a : 
sin {(2) F (fl + (le)*)'*"(t — wo} 
@ = Pe SP" (5.3) 


As further specialization of the above solution, consider the case of a flat plate by 
letting R — © and setting F = 0. Thus, for an infinite elastic circular plate subjected 
to a pulse concentrated at the origin (allowing a — 0 while P = 2a*p,), (5.1) reduces to 


Pp y\'/2 sin {(2) F(t _ tp 
(2) | tJ (rt) ——-=—_,-—_—_—+ & (5.4) 


wr, t tee 
2rph \D é 
which, with 4, = 0, the notation r, (Z/p)'’*t, and the use of [17, Table 8.2)] becomes 


6 p= SE) ee = 2)". Or (r)() — 
wr, t) = i ph? | 120 v’) 2 va | r i) \r, fr 5.5a) 
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si(r) = -| = = dx = —}tr + Si(r), (5.5b) 


Si(r) being the sine integral. 
In the remainder of this section, we confine attention to shallow elastic spherical 
shells subjected to loading of the types 


p(r, t) = —pH(a — re", (5.6a) 
pr, t) = 0 A(r)e'**, (5.6b) 


which vary harmonically in time. For such steady state solutions, it is more convenient 
to return to the viscoelastic solution (4.5) in the Laplace-Hankel transform-plane. 
Following Lee [6], we replace s by zw in all quantities except in p*(é, s) which is replaced 
by’ p*(£), and then take the inverse Hankel transform leading to the (real) steady state 
amplitude I(r) for the axial displacement. When p(r, ¢) is specified by (5.6a), through 


the process just described and with the aid of (4.14a), we obtain 
ee dre) J (ab) 
Wr) = —pa D' | “OR de, (5.7 
| Pol I ie + AS dé ‘) 


ith similar expressions for YV°f and df/dr, f(r) being the amplitude of F(r, ¢t) 
Unfortunately, the solution (5.7) does not admit a closed representation; and probably 
more direct approach, involving the use of the elastic solution of the 


rol th s Case 
problem [10] and the Laplace transform, will (from a practical point of view) prove 
fruit 


For the steady state solution due to an oscillating point load specified by (5.6b), 
the aid of (4.14b), W and V’*f are given by 


With 
. Pf” EJ o(ré) ia 
=__— - . - = ra hy 8: 
ti 2x D i fe" + ro] dé, (5.8a) 
V “f(r) : -*2 Wir). (5.8b) 


Che integral in (5.8a), with \§ as positive, admits a closed representation [17, Table 
8.2]. Hence, 
Wir) = - am ket(Xor) (5.9) 
(7) 7 Dr D2 ( 0 oO. 
which is in exact agreement with the results given in [10] for 
(R/y)? >1 or ww > (E/p)/R’ 


corresponding to Aj > 0); in (5.9), kei x (together with ber x, bei x, and ker x to be 
introduced presently) are the Kelvin functions. 

On the other hand, if in (5.8) 4 < 0, ie., for (R/y)? < 1 or w < (E/p)/R’, we 
introduce the quantity \ through 


No = \* exp (-:7) : (5.10) 


‘Here, this process is equivalent to taking the inverse Laplace transform. 
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[where, by (4.6a) with s replaced by iw, is defined by (3.4)] and also employ the relations 


[18] 
T ° ’ rr ul ae _ Z y & T 
kerz = ria exp ( i) aE exp ( *) a 
(5.1la) 
, i , {8 
+ iad exp ( *) —- aE exp (: *) | ' 
° T : ee a 7 
pie = 7{-a,: exp ( it)| + ard: exp ( it) —_— 
9] 11b) 


— 4 2 exp (i) | - iv.[2 ex (2) } 


whose argument z = x + iy is complex. Thus, for \§ < 0, substitution of (5.11) and (5.10) 
into (5.8) and the use of [15, Table 8.2] lead to expressions whose real parts’, i.« 


td , T yy, \ 3 , — 
W = Re - Dr | Ke +5 Yor) — 15 1.0.) | (5.12a) 
«cml. fe 
f =Re tr Dit (ME )[ wane) - 3 Yo(Aor) + is Jo(Aor) + 2 In “|| (5.12b) 


agree with the corresponding results in [8] for R/y < 1. Here, it may be of interest to 
note that with (5.10) and with the use of known relations of the type (see e.g. [18, p. 20]) 


Y,(i'°x) = [bei x + i ber x] + [—ker x + 1 kei x], 


Eq. (5.12) for the range R/y > 1 may be reduced to the form (5.9)*. 
6. Solutions for special viscoelastic materials. In this section, for simplicity’s sake 
we limit ourselves to incompressible media (v = 4), in which case (2.6) and (4.6) become 


E(s) = 3P.(s)P;"(s), 6.1) 


v, - 2) 1, os 
dels) = 1” E + a ‘ 


and deduce explicitly solutions for the viscoelastic Maxwell and Kelvin solids. 
(a) For the incompressible Maxwell solid with the aid of (2.7b), (6.1) reduces to 


E(s) = ——x E, 
slid (6.2) 
9 9 = E 


As tables for the inverse transforms of the functions involved here are available, we turn 
to (4.5) and write 


8It may be recalled that on account of the assumed form of the exponential time-dependence of 
the solution, the various quantities are in general complex, and their real parts give the desired results. 
In this connection, compare with [10, Eqs. (42) to (45)). 
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§ T ; A ‘ 
wr*(é, ss) = 3 2 sat ania Pp (é, $) 
Eh S ‘ 9p : ee 7 
; + Ek s +sr + R? 
. (6.3a) 
oa p*(é, 8) str’ 





oh s{(s + (27) ')? + (8° — (27)%)]’ 


B= (\(z) 0 + Chee (6.3b) 


rhen, for instantaneous pulse specified by (4.10), after recalling the Laplace transform 
of the 6-function [11, p. 27], « 


as well as its inverse transform [11, p. 323], and using [16, 
Table 5.2], we take the inverse (Laplace followed by Hankel) transforms of (6.3a) to 
obtain 


where 


1 


w(r, t) oh [ [éq*(£) J o(ré) ] = {1 — exp (- t— to r) 
] 1/2 
+ ‘OS Ue - 3) (t — 0) | (6.4) 





In a similar manner, 


’"F (r, t) A | léq*(é J 


pR Jo s o(r8)] (6.5) 
o(54) 
nice WE Wile 1\'? 
of —o x sin | (s" —_ i) (¢ = oj dé 
: peer 4r 
| (s 7 i) 
and dF /dr, as in See. 4, is obtained by integration. 
(b) For the incompressible Kelvin solid with the aid of (2.8b), (6.1) reads as 
E(s) = (1 + rs)E, 
6.6 
9p 8° + [E/(oR')\(1 + — 


No(8) = jy2 (1 + 7s) 


und again, to avoid the convolution integral in (4.7) and (4.9), we return to (4.5) and 
(4.8) and obtain 


w(r, t) = - i [Eq*(é) J o(ré)] exp | -2 ({ — «| 
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V'F(r, t) = = | [Eq*(E) J o(ré)] exp | -3 (t — 0) | 
: pk J, 2 


(i-(@Te-s) 
4 f : na - sin (sf - (2) | «WH ya 


bas | 


Tp 


Since the integrals in the solutions (6.4) and (6.5), as well as (6.7) and (6.8), converge 
rapidly, the numerical evaluation of the results is feasible. In addition, it may be noted 
that as tf ©, the second term in the integrand of (6.4) diminishes exponentially while 
the first remains finite. Hence, under instantaneous loading, the shell medium for the 
Maxwell solid will assume a permanent deformation; no such effect is present in the 
solution (6.7) for the Kelvin solid. That this observation is not unexpected becomes 
evident by merely recalling the absence and presence of a restoring force in the one- 
dimensional Maxwell and Kelvin models, respectively. On the other hand, for an oscil- 
latory load which itself (unlike an instantaneous load) supplies the restoring force, no 
permanent deformation takes place in the Maxwell medium. 

7. Shallow viscoelastic shell segments. The solutions for shallow viscoelastic 
spherical shell segments, in principle, may be obtained in a manner similar to those for 
unlimited shallow shells (Sec. 4); the chief difference, however, is the use of finite Hankel 
transform (in place of Hankel transform) together with some manipulations necessary 
to accommodate the edge boundary conditions of the shell segment. Although the 
method of solution (to be discussed presently) permits the treatment of shell segments 
with various boundary conditions, such as those considered in | 14, Sec. 1 1], for simplicity’s 
sake we confine attention to the case of shallow spherical shell segments (0 < r < ro) 
supported at r = rp by means of a ring which is restrained against rotation and non- 
resistant to axial force. Here, the regularity requirements at r = 0 are identical with 
(3.8), and the boundary conditions are given by 

ws, 1) = Niro, 2 = Me, t) = 0. (7.1) 

Before proceeding further, we recall that the finite Hankel transform of a (suitably 

restricted) function U’(r, s), as well as its inverse transform in the interval 0 < r < r, 


are defined, respectively, by [12, p. 83] 


U*(é; 3s) = | rJ (ré,) U(r, 8) dr (7 .2a) 
and 
9 J ort.) 
U'(r,s) = 3 UC: 8) ee (7.2b) 
To d [Ji (T&) 
where é,; are the roots of the transcendental equation 
J A7oe;) = 0 (7.2c) 


and in (7.2b), the summation is intended over all positive roots of (7.2c). 
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After introducing the new variables v and G 


— r: “ae 
— i-—_--—-————"y, , v% = V2 ATo , t ’ 7 ¢ 
’ u ie \ ), (7.3) 


G = VF — Gy, Go = VF (ro ; t), 
defined to satisfy 
our, , 8) = Vielro , t)} = Giro , ) = 0,7 (7.4) 


we return to the system of differential equations defined by the first of (3.2) and (3.5). 
Then, following the application of Laplace transform (with zero initial conditions) and 


with an appeal to the correspondence principle, we deduce 


saciranaith es | ” , r— 1, 
Dis Je’ + 5 (G' + G) p'(r,s) - ph a “e | . 


i 
hE(s) r— r; 
G’ = ‘ae ie ; v . 
R [ , ie | 


From the finite Hankel transform of (7.5), there eventually follows 


[-4 ' ar ‘ + , 
v*(t; , 8 g dv (s) | 133 Ao(s) J (oe, vs(s) 
” (7 .6a) 
) 
J (Tok) me. 8 
same dé ht 2 Gi(s) + e ~ ; +p , 
R Ds) &; Dis) 
hk(s) . | ‘ 
G*(£, , 8) 7 fe) + Aas) ]'§ —rod A (roe eis) 
. (7.6b) 
ha ae y*(E; , 8) 
= euros! G(s) + Eee 
R Dis) &, D(s) 
The two above equations may be put in the form 
v*(E; , 8) = voft(E; , 8) — Gofs(é; , 8) + vt(E , 8), (7 .7a) 
G*(E; ,s) = vigt(—; , 8s) — Gigt(é; , 8) + GE; , 8), (7.7b) 


where the functions f* , {* , g* , g , vt and G* are defined by direct comparison with 
(7.6). It may be noted here that, except for v* and G* , all functions involved in (7.7) 
are polynomials in the Laplace parameter s. 

In order to eliminate the functions vi and G/, we take the inverse Hankel transform 


of (7.7) and obtain 
v'(r,s) = vifi(r, s) — Gof(r, s) + v{(r, 8), (7.8) 
G(r, s) = vogi(r, s) — Gigi(r,s) + Gi(7, 8), (7.8b) 


and also by the second of (7.3), 


af 


r ok (r. 8) Gi n{l — gi(n, s)] dn + vi | ngi(n, 8) dn + | nG{(n, 8) dn. (7.8¢) 


or J0 J0 0 
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In view of the functional form of M, in (3.3), we next introduce the operator 
a= $=, (7.9) 


and with the aid of (7.3), by virtue of the last two of the boundary conditions (7.1) 
when applied to (7.8), determine vj and Gj as functions of ff , f} , gf , gi , vf , and GY, i.e., 


nt + [Afi(ro , 8) ][Afi(ro , 8)]" [ ngi an/ [ n(l — g) an} (7.10a) 


= [Afi(Tro 8) Aot(r ,8) — Afi(ro , 8) | nG{ an/ [ n(l — g5) an} 
70 v0 


\ 


and 


aro 


To —} ro 
Gs = -| f n(l — g2) in| {0 | ngi dn + [ nG{ in}. (7.10b) 
0 0 70 


Next, substituting (7.10) into (7.8) and taking the inverse Laplace transform followed 
by inverse Hankel transform, we obtain 


9 
o(r, t) = > DOL wiftles , 8); 
5: (7.11) 
e J o(ré;) 
— L“(Gifré, , 8); "wre; , 8); } —a 
ENGAGE: 5 8); + Le Wis, 8)5 a [J s(roé) |? 
|) . 
G(r, t) = = Do IL wigt(é. , 8); ¢] 


J ré;) 


L"[Gihgt(é; ,s): + L'(G*(é; ,s); t} — e 
[ Gls J | }} [Js (roé,) | 


which, together with (7.3) and (7.10), formally completes the solution sought. 


REFERENCES 
1, T, Alfrey, Non-homogeneous stresses in viscoelastic media, Quart. Appl. Math. 2, 113-19 (1944) 
2. E. H. Lee, Stress analysis in visco-elastic bodies, Quart. Appl. Math. 13, 183-90 (1955) 
. E. Sternberg, On transient thermal stress in linear viscoelasticity, to appear in Proc. 3rd U. S. Natl. 
Congr. Appl. Mechanics 1958 
4. M. A. Biot, Dynamics of viscoelastic anisotropic media, Proc. 2nd Midwest. Conf. on Solid Mechanics, 
pp. 94-108, 1955 
5. J. Mandel, Sur les vibrations des corps élastiques, Comptes rend. Acad. des Sci. 245, 2004-006 (1957); 
Also, Sur les vibrations des corps viscoélastiques 4 comportement linéaire, Comptes rend. Acad. des 
Sci. 245, 2176-178 (1957) 
6. E. H. Lee, Stress analysis in viscoelastic materials, J. Appl. Phys. 27, 665-72 (1956) 
. E. H. Lee, Viscoelastic stress analysis, Tech. Rept., ONR Contract Nonr-562(10), Brown University, 
July 1958; presented at the 1st Symposium on Naval Structural Mechanics 
. E. Reissner, On transverse vibrations of thin shallow elastic shells, Quart. Appl. Math. 13, 169-76 (1955) 
. K. Marguerre, Zur Theorie der gekriimmten Platte grosser Formanderung, Proc. 5th Intern. Congr. 
Appl. Mechanics, pp. 93-101, 1938 
10, E. Reissner, On azi-symmetrical vibrations of shallow spherical shells, Quart. Appl. Math. 13, 279-90 
(1955) 
11. R. V. Churchill, Operational mathematics, 2nd ed., McGraw-Hill Book Co., 1958 


wo 


“J 


© oo 





1960] SHALLOW VISCOELASTIC SPHERICAL SHELLS 121 


12 


13 


14, 


15 
16. 


. I. N. Sneddon, Fourier transforms, McGraw-Hill Book Co., 1951 

. J. P. Dalton, Symbolic operators, Witwatersrand University Press (Johannesburg), 1954 

E. Reissner, Stresses and small displacements of shallow spherical shells, J. Math. and Phys. 25, 

80-85, 279-300 (1945); 27, 240 (1948) 

B. Van Der Pol, and H. Bremmer, Operational calculus, Cambridge University Press, 1955 

Staff of the Bateman Manuscript Project, Tables of integral transforms, vol. 1, McGraw-Hill Book 

Co., 1954 

. Staff of the Bateman Manuscript Project, Tables of integral transforms, vol. 2, McGraw-Hill Book 
Co., 1954 

. W. Magnus, and F. Oberhettinger, Formulas and theorems for the special functions of mathematical 
physics, Chelsea Publishing Co., 1949 








122 BOOK REVIEWS [Vol. XVIII, No.2 
BOOK REVIEWS 


Introduction to functional analysis. By Angus E. Taylor. John Wiley & Sons, Inc., New 
York, Chapman «& Hall, Ltd., London, 1958. xvi + 423 pp. $12.50. 


The theory of linear operators and topological vector spaces grew from the study of the kind of 
linear differential equations and integral equations which occur in applied mathematics. Now that 
the pure mathematician has developed these theories, the applied mathematician is feeling the need to 
learn something of them. The publication four years ago of Bernard Friedman’s ‘Principles and Tech- 
niques of Applied Mathematics” clearly indicated the power of the new abstract methods in the analysis 
of practical problems. The book by Professor Taylor, under review, is the first textbook on functional 
analysis to appear in English (only the first part of the translation of the book by Kolmogorov and 
Fomin has so far appeared), and it is so well written that it will almost certainly become the standard 
text on the subject. 

The first two chapters ‘“‘The Abstract Approach to Linear Problems’ and “Topologies’’ provide 
the basis for the long third chapter on ‘Topological Linear Spaces’ which, despite its general title, 
is almost entirely confined to a discussion of normed spaces. Chapter 4, ‘““General Theorems on Linear 
is concerned with the study of continuous linear operators in normed linear spaces and 


Operators’, 
the principle of uniform boundedness, and 


contains basic results such as the ‘‘closed-graph theorem’’, 
the “projection theorem”’. The next two chapters,’’ Spectral Analysis of Linear Operators”’ and “Spectral 
Analysis in Linear Space’’, form, for the applied mathematician, the core of the book; the former deals 
with spectral theory in general Banach spaces while the latter is concerned with bounded operators 
in Hilbert space. The book ends with a chapter on integration and linear functionals, 

The whole book is exceedingly well written. In his preface, Professor Taylor writes: “This book 
is an introduction, not a treatise. It is meant to open doors for the student and to give him understanding 
and preparation which will help him to push on, if he wishes, to the new frontiers of modern mathematics, 
carrying with him a clearer realization of the structure of classical mathematics’’. This aim he has fully 
accomplished. New ideas are introduced carefully (and never without motivation) and many illus- 
trations of the applicability of the general theory are given. There are also problems for the reader 
to attempt. 

This is the most suitable book available for a graduate course on functional analysis and for in- 
dependent reading. An applied mathematician with little previous experience of abstract mathematics 
might find it difficult going. He had better begin by reading the book by Friedman or (if he can read 
Russian) Vulikh’s ‘“Vedenie v Funktshional’nii Analiz’’. But for the man who is well equipped and is 
prepared to make the effort, Taylor’s book is the most rewarding yet on functional analysis, though 


some pure mathematicians might think that already it is a little ‘‘old-fashioned”’. 
I. N. SNEDDON 


Handbook of automation, computation, and control. Edited by Eugene M. Grabbe, Simon 
famo, Dean E. Wooldridge. John Wiley & Sons, New York, Chapman «& Hall, Ltd., 
London, 1959. xxiii + 1096 pp. $17.50. 


The second volume of this unique handbook treats of: A. Computer Terminology, B. Digital Com- 
puter Programming, C. The Use of Digital Computers and Data Processors, D. Design of Digital 
Computers, E. Design and Application of analog Computers, F. Unusual Computer Systems. 

The applied mathematician will mainly be interested in Section B, in which John W. Carr III 
gives, in 270 pages, a masterly survey of programming and coding. About 150 pages are devoted to 
advanced topics, here presented for the first time in easily accessible form, principally the construction 
and use of automatic programming systems with examples drawn from existing U. 8. and Soviet com- 
pilers, interpreters and translators. The instruction logic of some common scientific and business com- 
puters, both large and small scale, is also given. The emphasis on existing computers and programming 
systems instead of hypothetical ones, often presented in books, is particularly welcome. 


WALTER F.. FREIBERGER 


(Continued on p. 130) 


123 


A ONE-DIMENSIONAL CONSOLIDATION PROBLEM WITH A 
MOVING BOUNDARY* 


BY 
R. E. GIBSON 
(Imperial College, London) 


1. Introduction. In a previous note [1] some simple one-dimensional diffusion 
problems of the sphere and infinite cylinder which involve a condition to be satisfied 
on a specified moving boundary were considered. The progress of consolidation in a 
clay layer increasing in thickness with time was discussed in a further paper [2]. In 
these examples it was supposed that the medium vanished at ¢ = 0, and therefore no 
initial conditions were prescribed; this restriction was made use of in the method of 
solution adopted. In many cases of practical consequence, however, a medium of finite 
extent exists initially through which the required quantity (e.g. temperature, concen- 
tration, water pressure) has a given distribution, and in the present note we shall be 
concerned with a problem of this type, taken from the field of soil mechanics. 

2. Consolidation of a clay layer. Consider a layer of homogeneous saturated clay 
bounded by the planes x = 0 and x = h(t). The excess water pressure u in the pores of 
the clay [2] is governed by the equation 


du _ du dh (1) 


. ax” ~—s tex Y dt ’ 


where ¢ and y are respectively the coefficient of consolidation and the submerged density 
of the clay. We seek a solution of (1) subject to the conditions 


oe (0,1) =0, t>O0, (2) 
ufh(t), t]} = 0, t>0, (3) 

and 
u(x, 0) = f(x), 0 <2 < h(0), (4) 


where A(t) and f(x) are given functions. 

A direct approach, using for example a Fourier type integral transform, would 
appear to carry us no further than can the following elementary considerations. A 
solution of (1) appropriate to the conditions (2) and (4) is clearly 


h(0) 


wa, t) = y[h(t) — h(O)] + (4xet)'”” [ f{(E) K(x, t; &) dé, (5) 


v0 


c—e (x + £)’ 


*Received February 19, 1959; revised manuscript received June 12, 1959. 


where the kernel 
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To this we must add a solution to (1) which will enable the requirement (3) to be met 
while not disturbing conditions (2) and (4). A function possessing the necessary proper- 
ties is 


v(x, t) = —(4mct)™'”” [ g(t) K(x, t; &) dé, (6) 


7A(O) 


since 


v(x, O+-) — g(x), h(0) <2 < @ 


= Q, 0<2z < A(0), 
where the function g is the arbitrary element we need to satisfy (3). It is, of course, 


assumed that the integral (6) exists for ¢ > 0; this must be verified a posteriori. 
The problem therefore reduces from (3), (5) and (6) to the following integral equation 


=) 


w{h(t), t] = (4ect)"”? | 


g(E)A[h(t), t; €] dé (7) 
7h(0) 
for g(é). For reasons which will appear presently it will prove convenient to write 
9=n + 9 


and replace (7) by two integral equations, namely 
y{h(t) — h(O)] = (4xct)7'”” | gi(&)K(h, t; ) dé (8) 
#7A(0) 


and 


=] 


[- f(t)K(h, t;&) dé = | 


A(O 


g(&)K(h, t; &) dé. (9) 


3. Reduction of the first integral equation. It is unlikely that these equations can 
be solved for an arbitrary A(t), and we shall be concerned here with the following simple 
relation 


h(t) = a+ mt 


where m, the velocity of this upper boundary, is constant, and a > 0. 
Following Parodi [3] we take the Laplace transform of (8) which becomes 


2 o 2 
(4xc)'?ymp~* = / gi(é) / i’? exp — E + mite a") dt dé 
a 0 ; 


+ | g;(&) [ t’”? exp — E + (mt tate dt dé. 


With the following changes of variable 
ec" rT¢, 


“ww 
a (p + m) 


& 
| 
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and 
G(r) = g(r + a) 


and proceeding formally, we find that 


2\—<- 


; 8 m , 8 m 2s m 1/2 ( 2 m ) . 
7* —-= G* 5 { jn — oe ; —la —_ rym —_— «<< 0 
(; 7 *) (; ; =) = | (a ”) | ne 4c} ’ On 


7*(d) = fg e-™'G(r) dr. 


We now put 


in (10), which becomes 
2 


2 -2 
G*(z) + re + m) exp | -20(¢ + ") = 2yme( + ml (c +. m) c- mn] », (I) 


and in this way the integral equation is reduced to a linear difference equation. 
4. Solution of the difference equation. The difference equation may be treated by 
considering first the solution of (11) with the right hand side zero: 


o(z) + 6(: + =) exp | —2a(s + | = 0; 


4 


using elementary methods [4] we find 


¢(z) = exp | a1 + e) + oe (12) 


m J 





as a special solution. To obtain the required solution put 
G*(z) = o(z) Wz) (13) 
and it follows from (11) that 


(2 + ™) — #(s) = —7 = (2. + m).-(s + m) exp {—[ a + ) + ire, |} 


so that, formally at least, a solution is 


We) = 7 > {  - E + (n + 1) my)" 
a=0 ( "| c (14) 


-exp {-[ (2 + nm +n+ «) + in(n + #) |. 


Using equations (12), (13) and (14) the following expression is obtained 


og =r Eco] - [+42] ] 
“ (15) 
“exp {—an| 2 + (n + 1) m | 
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5. Solution of the integral equation. The required function G(r) of Section 3 can 
now be found by inverting (15), and we use the result 


if : ; , 0, O<r< 
| ee + 8) de = J S7SP 


271 
(r-—Be "7" T]Bs 
By introducing the unit function 
U(r) = 1, r>O 
= 0, 7 <0 


the solution can be expressed in the succinct form 


G(r) = ¥ >> (—1)'(7 - 2an)\exp | — * 7 = Zan) 


. m . am 
— exp | —(n + 1) — (tr — 2an) b (7 — Zan) exp | —— n(n + ] ; 
‘ c } Cc 
| 


: ‘ f WIR ,. 
gil(é) = ¥ > (-1) [fé — (2n + l)a}§ exp |= (§ — 2na — «| 
; ; c 


so that, finally, 


(16) 
m \ . am 
— exp| m+ (¢ — na — a) |Sve — @n + 1a} exp| —* “n(n + | | 
c ) c 
The water pressure u’(x, ¢) corresponding to an initial pressure f(x) = 0 in the layer 


0 < x < a, the pressures being generated solely by the increasing weight of super- 


incumbent material, is given from (5) and (6) by 
u(x, t) = ymt — (4qet)'”” [ gi(é)A(x, t; &) dé. 
Rewriting (16) in the form 
g(t) = > T(é, n) UE — (2n + Ia], 
we have 


u’ = ymt — (4acl)™' [| T(é, O)K (x, t; &) dé + | T(é, 1) K(x, t;£) dé + | 


5a (17) 
= ymt — (4nct)'” >> | T(é, n) A(x, t; &) dé 
n=0 J (2n l)a 
and the integration may be accomplished by introducing the variable 
¢=£& — (2n+ la 
in (17); we find finally that 
u’ = ymt — y(et) ~ exp (2 :) > (—1)"[F(az, t,n) + F(—a, t, n)] 
cl) 
(18) 


[ am a” =) 
“exp j— n(n + 1) + —-(2n + 1)* |? 
c tet 
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where 
F(x, t,n) = [x(a2) — x(a,)], 
x(a) = ae*’ erfea, 
a, = (4ct)~'”*[2nmt + (2n + loa — 2], 
and 
a, = (4ct)"'”*[2(n + 1)mt + (2n + la — a]. 
We note here an earlier result [2] for the case a = 0: 


> 


. o © — 
u'(x, t) = ymt — y(ret)'’* exp (-=) [ € tanh cosh x exp (=) dé. 


The more general solution (18) with a = 0 can easily be shown to be equivalent to this 
if the term sech (mé/2c) in the above integrand is expanded in exponentials and the 





series integrated term by term. 

6. Reduction of the second integral equation. In the above treatment of the first 
integral equation (8) we found it possible to reduce the problem to one of determining 
the solution of a difference equation in the complex plane. The solution of the second 
integral equation (9) can be found in the same way, but it is instructive to employ a 
direct method by which the problem is reduced to the solution of a difference equation 
in the real plane. 

We rewrite equation (9) in the equivalent form 


[ f() exp [—(mt + a — £)?/4et] dt = / g(t)K(a + mt, t; ) dé (19) 


v—a 


where {(£) is a specified even function. 
We now impose the condition 


g(t) = 0, -a<t<a (20) 


and the right hand side of (19) then becomes: 


mt\ f° , . __ (mdr : md : ? 
exp (- te | Ee + a) exp (™) + g{XA — a) exp ( Qe )| exp ( -) dx, 


where we have put £ = \ + ain the first term, and = \ — a in the second term of the 
integrand. The left hand side of (19) is now modified by taking 


f(é) = 0, || >a, (21) 


and then 


[ f(z) exp [—(mt + a — £)’/4et] dé = | f(é) exp [—(mt + a — £)*/4ct] dé 


= [ f(\ — a) exp [—(mt + d)’/4ct] dd 
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where \ = a — &. The equation (19) now takes the form 


[ Ee + a) exp (=) + {g(A — a) — f(A — a)} exp (-™)| exp (-.) Ii = 0 


for all values of ¢ > 0. It is natural to assume that the integrand vanishes and this leads 
to the following difference equation 


go(\ + a) exp (md/c) + gA — a) = f(A — a). 


The condition (20) enables g, to be determined uniquely. We omit the details of the 
solution which is 

g(t) = (—1)"*'f(— — 2na) exp [—nm(~ — na)/c); (Qn —1) <= < (Qn + 1). (22) 
From (5), (6) and (22) the water pressure u’’(z, ¢) in the layer 0 < x < a + mt, when 
the initial pressure is f(r) in 0 < x < a and when the accreting material is supposed 


weightless (y = 0), is 


u’"(x, t) = (4ect)”*”* | f(é) exp [—(x — €)*/4et] dé 
aa (23) 
— (4ret)’”” 7. (—1)""" | f(é)K (x, t; & + 2na) exp [—nm(E + na)/e} dé 
The complete solution to the problem of Section 2 is now found from the sum of (18) 
and (23) 


u=uy’ + y’’. 


It can easily be verified that equations (1)—(4) are satisfied and that the series involved 
in the solution converge for t > 0 and f bounded. 

7. Some alternative boundary conditions. When other boundary conditions are 
imposed the use of the above method depends upon our being able to choose a solution 
at the outset which satisfies all conditions except that on the moving boundary; in simple 
cases a knowledge of the appropriate Green’s functions will be sufficient. For example, 
to satisfy (1) with initial condition (4) and boundary conditions 


u[h(t), t] = (2) (24) 
u(O, t) = ¢(b) (25) 


we should take as a tentative solution 


u = y[h(t) — h(O)] — 2(4ae)'” [ Pit NO) Ko) exp | -is| dr 
Jo 7 , 7 (26) 
fo) A(O0) 
+ [ g@Q@, 59+ [  1OQ@, t;8 ae, 
~ ACO) 70 
where 


Q(x, t;&) = (4nct)'’* {exp [—(a — £)/4et] — exp [—(x + £)*/4ct]} 


S 
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our choice being dictated by the nature of the initial and boundary conditions. The 
expression (26) satisfies (1), (4) and (25) and the function g is to be found from the 
condition (24); this requires us to solve the linear equation 

' fyh(r) — yh(O0) — $(7)] | _ hb) 5| 


Hi) — 1h) — MO) + Hoey" | 7 Heap | 
Y(t y| h(0) ] h(t)(4xc) - i 7)? ex 4c(t — 7) 


—f[ samo, ceae= [ «QMO, t; 8148, 


70 


and with /(t) of the form (10) the solution may be found along the lines adopted above. 
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Einfiihrung in Theorie und Anwendung der Laplace-Transformation. By Gustav Doetsch. 
Birkhauser Verlag, Basel and Stuttgart, 1958. 301 pp. $9.15. 


Professor Doetsch’s three volume Handbuch der Laplace-Transformation has established itself in 
the last ten years as the standard treatise on the Laplace transformation and its applications. The 
volume under review aims at being a short introduction to the subject. The book is divided into 28 
sections with chapter divisions. In the first eleven sections the elementary properties of the Laplace 
transform are developed. Sections 12-15 describe the application to ordinary differential equations 
with constant coefficients, systems of ordinary differential equations and difference equations, while 
§16 is concerned with the behaviour of the Laplace transform at infinity. The central part of the book 
(§§17-25) is devoted to a discussion of the inversion formula and related topics including Parseval’s 
theorem and the asymptotic behaviour of the image function and the original function. The final three 
sections are concerned with applications: §26 considers the solution of ordinary differential equations 
with polynomial coefficients (in particular, Bessel’s equation), §27 deals with partial differential equa- 
tions in two independent variables (the diffusion equation and the telegraphy equation), and the last 
section is given over to a discussion of the solution of integral equations of the ‘‘Faltung”’ type 

The book is clearly and attractively written and anyone wishing a sound introduction to the theory 
of the Laplace transform is well advised to read it. For all that, an applied mathematician or engineer 
beginning the subject would probably find Churchill’s ‘Operational Mathematics” more to his taste. 


I. N. SNEDDON 


Theory of differential equations. By A. R. Forsyth. Dover Publications, Inc., New York, 
1960. Vol. 1, xiii + 340 pp.; Vol. 2, xi + 344 pp.; Vol. 3, x + 391 pp.; Vol. 4, xvi + 
534 pp.; Vol. 5, xx + 478 pp.; Vol. 6, xiii + 596 pp. $15.00. 

This unabridged republication of the well-known work consists of three volumes each of which 
combines two volumes of the original edition. 


Theoretical elasticity. By Carl E. Pearson. Harvard University Press, Cambridge, 1959. 
218 pp. $6.00. 


This book evidently grew out of the author’s lectures given at Harvard University, the purpose 
of its publication being “‘ in part to discuss modern methods of elasticity theory in a form which 
does not require extensive mathematical experience and in part to provide a compact and convenient 
summary of such methods.” The selection of topics discussed deviates from the usual coverage in books 


on elasticity. Thus, for example, plane problems and Saint-Venant’s torsion problem are excluded. 
On the other hand, a few of the topics dealt with in the latter part of the book (notably parts of chapters 
6, 7, and 10), may not be found in other books dealing exclusively with elasticity. No attempt is made 
to supply adequate and appropriate references to the brief treatment of the various topics; a few refer- 
ences are cited as footnotes. The author is unusually consistent in his omission of references, however, 


since for example, a good deal of the content of the last chapter is based on the author’s own excellent 


paper {Quart. Appl. Math. 14, 133-144, 1957] to which he does not refer. 


Chapters 1 and 2 (45 pages) are concerned with a review of the vector algebra, vect ileulus, 
index notation, and Cartesian tensors. Chapters 3, 4, 5 (63 pages) are respectively devoted to the analysis 
of stress, analysis of strain, and derivation of the basic equations together with some of the funda- 
mental theorems in the linear theory of elasticity. The next two chapters (6 and 7) are entitled ‘General 
Solutions’ and ‘‘Variational Methods,” respectively. The former, following a brief review of potential 
theory, contains an account of various integral representations (such as Betti’s formula for dilitation) 


(Continue d on Pp. 154) 
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PLATE DESIGN FOR MINIMUM WEIGHT* 


By 
R. T. SHIELD 


Brown University 


Summary. A discussion of the equations determining the minimum weight design 
of transversely loaded sandwich plates is given. By means of an inverse method, the 
solution for an elliptical plate clamped around its edge is derived. 

1. Introduction. Apart from a paper by Prager [1|f, previous work [2-5] on the 
minimum weight or minimum volume design of plates has been confined to circular 
plates. In the present work, no assumption of circular symmetry is made. The basic 
equations are formulated and the four different types of solution are discussed. An 
inverse method is developed to obtain the minimum volume design for plates with 
built-in edge conditions. The method is used to obtain the design for an elliptical plate. 

2. Basic equations. It is required to design a plate of given shape and with prescribed 
conditions of simple or built-in support along its edge. The design is to be such that 
the plate supports by bending stresses a given distribution of pressure and in addition 
the consumption of material is to be kept to a minimum. Here it is assumed that the 
plate is a sandwich plate consisting of a core of a given constant thickness H and identical 
face sheets of (variable) thickness h, where h «< H. The core carries shear force only 
and a bending moment across a section of the plate is supplied by direct stresses in the 
face sheets. The material of the face sheets is assumed to be elastic-perfectly plastic 
and to obey the Tresea yield condition of constant maximum shearing stress during 
plastic flow. The principal bending moments 1/, , /, at a point of the plate are then 
restricted to lie within or on the hexagon of Fig. 1 [6]. The maximum bending moment 
M,, is equal to o)fh, where op is the yield stress of the face sheet material. Plastic flow 
ean occur only for states of stress represented by a point on the yield hexagon, and 
during plastic flow the normality condition applies. This states that the vector with 
components proportional to the principal plastic curvature rates x, , k2 is normal to 
the side of the hexagon for stress points on a side, and lies in the fan bounded by normals 
to the adjacent sides for stress states represented by a corner of the hexagon. 

Because of the normality condition, the rate D of dissipation of energy due to plastic 
action per unit area of the middle surface of the plate, given by 


D = Mn, + Max, , (1) 
is uniquely determined by the curvature rates x, , ko . It has been shown [3, 7, 8] that 
minimum volume of the face sheet material will be involved if the plate is designed to 
be just at the point of collapse under the given loads in a deformation mode such that 

D/h = constant (2) 


‘Received June 23, 1959. The results presented in this paper were obtained in the course of research 
sponsored by the Office of Ordnance Research, Department of the Army, under Contract No. DA-19- 
020-OR D-4795. 

tNumbers in brackets refer to the list of references at the end of the paper. 
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Fic. 1. Yield condition. 


over the plate if the weight of the plate is neglected. The result is not restricted to plates 
in bending but also applies to shells [9] and weight can be taken into account. The term 
collapse is used here, as in limit analysis [10], to denote conditions under which plastic 
flow can occur under constant loads, changes in geometry being neglected. It should 
be noted that the dissipation rate D is directly proportional to the face sheet thickness 
h so that condition (2) does not involve h, and is a condition on the plastic curvature 
rates x, , Kk, only. However the form of condition (2) does depend upon the position of 
the stress point on the yield hexagon. 

During collapse, the elastic strains are constant [10] so that the plastic curvature 
rates are determined directly from the transverse deflection rate w of the middle surface 


of the plate by the equations 


ow , 1 dw | 

1 oe 5 ate | 
ad 1 e2 fe 
4 ) 3) 
Ow , 1 dw | 

me =-~r+—--.| 
OS P2 OS; J 


In these equations, s, and s; denote distance along the (orthogonal) lines of principal 
curvature rate which coincide with the lines of principal bending moment. The quantities 
p, and p, denote, with due regard to sign, the radii of curvature of the lines of principal 


bending moment, Fig. 2, and we have 


= 0¢ , dy 
i = 8 a a 9 (4) 
Os Os : 


where ¢ is the angle between the first principal moment direction and the x-axis. The 


radii p; , p2 satisfy 
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! 
a 
Fic. 2. Principal directions. 
and also 
a 1 @ a 1 @ . 
——-1#__-%__142. 6) 
08, O82 p, O8, OS» OS, P2 O82 


Because the s; , s, directions are principal directions the curvature rate x. is zero, so 
that 

Ow 1 dw a°w 1 dw (7) 

-=- - —— —_-——- d 





OS, O82 p; Os,’ OS, O8; Po O8, 


Equations (3)—(7) can be used to show that 


\. 
<1 = . (ke — &), (8) 
Os p; 
IK» ] 
2 = (Kk, — Ko). (9) 
Os, Po 


As these equations are obtained by eliminating w between Eq. (3), they are the com- 
patibility conditions on xk; , ke. 

The bending moments M, , , and shear forces Q, , Q. per unit length acting on an 
element of the plate with faces normal to the principal directions are shown in Fig. 3. 





Fic. 3. Positive senses of shear forces and moments. 
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For definiteness we assume that the plate is horizontal and that the pressure p acts on 
the upper surface of the plate. Positive values of 17, , M, stress the lower face of the 
plate in tension and the shear forces are positive in the senses shown in the figure. To 
correspond, the deflection rate w is positive in the upwards direction. For equilibrium 


aT AO), 

3, + <2 + : Q, + Q2 = Pp, (10) 

OSs, O82 Pez Pi 

and 

aM | 

SS te, + 9 +G, = 9; 

ds, P2 | (11) 
| 

7] : 

4 La, - uM) +0, = 0.) 

O82 Pi , 


{eferred to the rectangular cartesian coordinates (x, y), the equations corresponding 
to Eqs. (3), (10), (11) are 


Ow Ow ow 
kK a hy aes Kry i oe (12) 
ON OY Ox OY 
00), Ow), . 
: a ae, (13) 
OX OY 
OM, aM, P 
— — + @, = 0, 
OX 01 
Y ( | }) 
aM, . aM 
: + — + YY, = 0. 
OX OY 
Equations (13) and (14) can be combined to give 
aM, a’°M,, . aM, x 
2 ae: P - eo 4, .2 —~?p- (15) 
Ox Ox OY oy 


3. The four types of solution. The condition (2) for minimum volume imposes a 
restriction on the deformation w. The form of the restriction depends upon the position 
of the stress point on the yield hexagon of Fig. 1. There are only four essentially different 
types of plastic flow which can occur. For example plastic flow with moments represented 
by corner A of Fig. 1 differs in character from plastic states at D by only a change in sign, 
and stress points on BC become stress points on /F if the suffices 1, 2 are interchanged. 
As typical of the four types of plastic flow, we shall consider in turn moment states 
represented by corners A and C and points on the sides AB and BC, pot including the 
end points. For convenience in the algebra, the constant in condition (2) will be taken 
throughout to be ko,H, where k is a positive constant. It will be seen that regimes AB 
and BC lead to restricted forms for the deflection rate and it is unlikely that these 
regimes will play a large part in the solution of particular problems. 

Regime A. For the corner A, M, M, = M, = oo ffh and the curvature rates «x, 
and «x, are positive. The condition (2) then requires x, + x. = k, or 

Ow aw 
— + —, = k, (16) 
Ox” oy ; 
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Equilibrium requires 


aM, , #M, 


Ox" + ay” ae (17) 
Regime C. For the corner C, M, is zero and M, = — M, . The curvature rates must 
satisfy the inequalities 
—K, > x, >0, (18) 
and the constancy of D/h requires 
x, = —k. (19) 
rom Eqs. (8) and (19) it follows that 
(xo + k) = 0. (20) 
pi 
\s ko is positive and & is by definition a positive constant, this equation requires p) = © 


and the lines of principal moment in the s,-direction are straight. Thus the orthogonal 
net of principal moment trajectories consists of families of straight lines and parallel 
curves. 

[t is convenient at this point to introduce new coordinates, Fig. 4. One of the curved 





Fic. 4. Coordinate system for regime C. 


principal moment lines is chosen as a base curve and distance from this curve along 
the straight lines is denoted by «. The inclination ¢ and the distance @ are used as co- 


ordinates, and we have 
ds, da, ds, : Pe de, pP2 = Pp + a, (21) 
where p p(v) is the radius of curvature of the base curve a = 0. Equation (19) is 
simply 
a* p 
<——; = —k (22) 


Oa 
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and therefore 
w = —tka + af(y) + g(¢). (23) 


The functions f and g are not completely arbitrary but must be chosen to satisfy Eqs. 
(7). The final result is 


w= —tka® + af(y) + | p(y) f(y) de, (24) 


d 


where the prime denotes differentiation with respect to ¢. In addition, inequalities (18) 


must be satisfied. 
It is easily shown that the equilibrium equations and the yield conditions MV, = —M,, 


M, = 0 give 


IM, M, 
Q, = aM ab. ‘ 0. = 6, (25) 
Oa Po 
D(¢) : tf s 
i «tg 2 « | ppda ~ | p2p da, (26) 
Pe d P2 


where C and D are functions of ¢ only. 
In the particular case when both families of principal moment lines are straight, it 
is convenient to take the z, y-axes in the s, , s,-directions. We then have 


aw ow 

Pe —k, ao = U, (27) 

Ox Ox OY ‘ 
so that 

w= —tkrx + ax + gy). (28) 


Here a is a constant and g(y) must be such that 


k>qg'(y) = 0 (29) 


in order to satisfy inequalities (18). 
Regime AB. For moment states represented by points on the side AB of Fig. 1, 


M, = M, and M, > M, > O. The curvature rate «x, is zero and condition (2) then 
requires x. = k. As both x, and x, are constant, Eqs. (8) and (9) require p, = pz oo 
and both families of principal moment lines are straight. If the xz, y-axes are taken in 
the s, , s.-directions, we have 

Ow ow Ow 

— = - 0, “2 (30) 

Ox OX OY OY 
and therefore 

w= thy +ar+ byt+e (31) 


where a, 6, ¢ are constants. 
The moment M,, is zero, M, = M, and M, > AM, > 0, and the equilibrium equation 


(15) requires 


} 
| 


aM, , aM 


— = +f : = —®f. (32) 
i Oy / : 
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Regime BC. In this case x, = — xk. = — kand again both families of principal moment. 


lines are straight. With the x, y-axes parallel to the s, , s,-directions, the deflection rate 


is given by 
w= —hk(a’ — y’) +axrt+ by +e. (33) 
The moments satisfy 


M,, = 0, M,-—-M,=M,, M, > M,> 0, (34) 


“ zy 
and for equilibrium 


aM, , aM, 

ox - ~ —p. (35) 
4. Simply supported plate. When the plate is simply supported at its edge, the 
deflection rate w and the normal bending moment M, must be zero at the boundary 
of the plate. For a circular plate, the minimum volume design is obtained by the use of 
regime A only. This is also true for a certain class of shapes. If the plate shape is such 
that the deflection rate w determined from Eq. (16) and the condition w = 0 at the edge 
gives non-negative curvature rates x, and x, , then the deflection rate satisfying D/h = 
constant has been found. The thickness distribution for any given distribution of pressure 

is obtained by solving Eq. (17) for M, subject to M, = 0 at the edge. 

lor example, for the plate bounded by the ellipse 

ae 


;=1, (36) 


ce: 
the solution of (16) satisfying w = 0 at the edge is 


hath? (x yy ) . 
= 2(a” + 6°) (= b° I), (37) 


and x; , k2 are non-negative. For constant pressure over the plate the thickness distri- 
bution is given by 
pa*b* fo 
ae re (1-5 -§) (38) 

The case of the simply supported plate has been discussed previously by Prager [1]. 

5. Built-in plate. For the plate with a built-in or clamped edge, the deflection rate 
w and its normal derivative dw/dn must be zero at the edge. The edge is therefore a line 
of (zero) principal curvature rate. The solution for the circular plate involves regime A 
in a central circular region and regime C in the remaining annular region [5]. Here an 
inverse method is developed for obtaining solutions similar in character to that for the 
circular plate. 

It is assumed that in an inner region, which is bounded by the curve y, regime A 
applies, and in the region between y and the edge I of the plate, regime C applies, see 
Fig. 5. The curve y is as yet undetermined. The lines of principal moment in the region 
between y and I are straight lines normal] to T and curves parallel to I’, since the boundary 
is a line of principal moment. The trajectories are indicated by the dotted lines in Fig. 
5. The coordinate system of Fig. 4 will be used, with the base curve a = 0 taken to be 
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Fic, 5. Plate with built-in edge. 


the edge T of the plane. The direction of a increasing is taken in the direction of the 
outward normal to I, so that @ is negative in the region between y and I. 

The deflection rate we in the annular region has the form (23) and in order to have 
w = Oand dw/da 0 on T (where a = 0) we must have 


We = —htka’. (39) 


The curve vy and the deflection rate w, inside y are to be determined from Eq. (16) 
which holds inside y and the continuity of w and dw/dn across y. Assuming that the 
solution exists (with non-negative x, and x inside y), the determination of y may prove 
to be difficult, and for this reason an inverse method of solution is developed here. 

From Eq. (39), we have 

OW, 
- = —ka = (—2kw,)’. (40) 
da 
Also since we depends on @ only, dw-/da is equal to the magnitude of the gradient of 
w- . It follows that the deflection rate (39) is such that 


Jw, \° Aw, \" 
Qkwe + (— + (= ) = 0. (41) 
Ox OY 


This equation holds in the region between y and [ and in particular Eq. (41) holds on 
the curve y. Thus the deflection rate w, satisfying Eq. (16) inside y must satisfy (41) 
on y because of the continuity of w and its first derivatives across 7, that is 


; dwal’ , (dw, \" 
Qkw, + (2 s) _ (< +) = (0 on y. (42) 


O21 ; OY 


The inverse method is therefore to choose a function w, which satisfies Eq. (16). 
The curve y is then determined as the locus of points satisfying Eqs. (42). In order to 
determine the curve I which defines the boundary of the plate, we again use the con- 
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tinuity of w and its first derivatives. From (39), the (variable) value of a on y is de- 


termined by 
. 1/2 
ae (- ‘) (43) 
k Y 


where the suffix y means that the value on y is to be taken. Also in the region between 
y and I, the inclination ¢ of the straight principal moment trajectories is given by 


dw./dy 
ms =———— 44 
° OW /dx ’ (44) 
since these lines are normal to the curves w = constant. Thus the value of ¢ on y can be 
determined from the values of dw,/dx and dw,/dy on y¥, 
Ow,/01 ; 
tan g, = (2" L0H) : (45) 
OwW,/Ox/ , 


Knowing the values of a and ¢ on y, the curve I’ can be found. The coordinates of a 
point on y will be denoted by (&, 7) and the point on I’ obtained by dropping a perpen- 
dicular from (£, ») to T will be denoted by (X, Y). The points (£, ») and (X, Y) then are 


distant — a, apart and the line joining them is a principal moment line. It follows that 
X =t- a, cos¢,,) 

¢ / Py ‘ (46) 
Y = 7 — Qa, sin Ly J 


With Eqs. (42), (43) and (45), Eqs. (46) can be written as 


= ] (2a) 
=f = | 
- a oy 


1 2a) 
Y = n+ i(@ ? 


and these equations are the parametric equations of the curve I, remembering that 


(47) 


| 


the point (é, 7) lies on y. 
It should be remembered that «x, and x, must be non-negative inside y and the in- 
equalities (18) must be satisfied in the annular region. With w given by Eq. (39), 


ka ; 
ie —&k, aa (48) 
p> 
and since a is negative the inequalities (18) will be satisfied provided that between 


y and IT, 
w= i. (49) 


It is only necessary to check this inequality on y as the left hand side attains its maximum 
on y. Remembering the definition of p, [Eq. (21)], inequality (49) can be interpreted as 
requiring that the radius of curvature of T at (X, Y) be not less than twice the distance 
— a, between the points (£, ») on y, (X, Y) on I. 

The thickness distribution, or equivalently the values of M, , throughout the plate 
can be determined for any given distribution of pressure when the curves y and T are 
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known. At the junction y of the regimes A and C, equilibrium requires the bending 
moments M, , M,, to be continuous and the shear force Q,, must be continuous. In order 
to make the bending moments continuous, M, must be zero on y. This condition and 
Eq. (17) then determine M, in the region bounded by y. The value of the shear force 
@, on y can then be found. The values of J/, in the annular region are then given by 
Eq. (26) where the functions C(¢g) and D(g) are determined by the condition /@, = 0 
on y and the known value of the shear force on y. 

We remark that if the pressure distribution p(x, y) is non-zero at points inside y, 
then the thickness h is non-zero everywhere in the plate except for points on the curve 
y. However, if p is identically zero inside y, then h is zero in this region, as h satisfies 
Laplace’s equation inside y and is zero on y. In effect, the load produced by pressure at a 
point in the annular region is carried by a “cantilever beam” from the point to the 
nearest point of the edge of the plate. Loads in the inner region are carried by the whole 
boundary of the plate. This property of the minimum volume design for a built-in plate 
has been discussed previously for the case of the circular plate [5]. 

6. Elliptical plate with built-in edge support. If it is assumed that inside y the 
deflection rate w is given by 


bol 


kr? 2 x 2) 
Wa = —$ a {1 -=-4 (50) 


(+9) ' 


. 2 ? 
7 q 


where 7 and g are constants, the solution for the elliptical plate is obtained The de- 
formation rate (50) satisfies Eq. (16) and x, , x, are positive. The curve y is determined 
by substituting (50) into Eq. (42) and the result is the ellipse 


where the semi-axes a and 6 are given by 
27,2 2 2/2 2 
» r(r+q) p2 qr +q) (59) 
= . _ — ‘ ) = = = ° OZ 
r+ 2¢° 2ar+q ) 


The values of a and ¢ on y can be obtained from (50) and (43), (45). The substitution of 
(50) into Eqs. (47) gives the coordinates (X, Y) of a point on T in terms of the coordi- 
nates (&, ») of the point of intersection of the curve y with the normal to T at (X, Y), 


E(r? + 2q’) y _ n(2r° <4. q). (53) 


wae +) ’ (r° + q’) 


As the coordinates (£, 7) satisfy (51), the equation for the boundary I of the plate can 


Ss 


be written 


7 } 2 7 = 
I . A? + RB’ = ly (54) 
where 
: r(r? + 2q°) 2 q (2r° + q) an 
A —__ —- B’ = te 5E 
fre? ae (55) 
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Equations (55) can be solved for r’ and q’ to give 
3r° = 2A? — B’ + (A* — A*B’ + B)'”, 
3q° = 2B’ — A*® + (A* — A°*B’ + B)™. 
The deflection rate we at the point (x, y) in the annular region is given by 
We = — ka’ = —hk[(X — x)? + (Y — y)’], (57) 


where the normal to T at (X, Y) passes through (2, y), Fig. 5. Because the line joining 
(x, y) to (X, Y) is normal to I, 


(56) 


Y-y_ YA 
X-—-x XB 


For given (x, y), X and Y can be determined from Eqs. (54) and (58) and substitution 

in (57) then yields the value of we at (x, y). It can be shown that inequality (49) is 

satisfied in the annular region, verifying that regime C applies in the annular region. 
In Fig. 6 the values of a/A and b/A are plotted against B/A for the range 0 < B/A <1. 


(58) 
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Fic. 6. Semi-axes a and b of curve y for elliptical plate with semi-axes A and B. 


When the plate is circular (A = B), the curve ¥ is a circle of radius two-thirds of the 
radius of the plate, as found previously [5]. For a plate which is long compared to its 
width (B/A — 0), the ratio a/A tends to unity and b/B tends to one-half. This can be 
compared with the minimum volume design for a sandwich beam, built-in at the ends, 
for which the sections of zero thickness, which correspond to the curve , are distant 
one-fourth the length of the beam from each end. 

Thickness distributions for specific pressure distributions may now be obtained as 
described in the previous section, and for illustration we consider the case of constant 
pressure p. The notation My, and Moc will be used to denote the value of My) = o.Hh 
in the inner and annular regions respectively. As M, is zero on y and M, satisfies Eq. 
(17) inside y it follows that the value of M, at points inside y is given by 


foe or a ¥). 
Moa ap (a° + b*) ] a’ b? (59) 
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In the annular region, W/, is given by Eq. (26) and as p is constant and J/ 0 on ¥ 
it follows that 
| , 2 I I 
Moc = = P(p2 — 13s) + Diy)\— — i (60) 
6 Pe2 P3 


where p,(v) is the value of p, on y. The value of D(g) is to be determined by the condition 
that the shear force Q, be continuous across y. It is convenient in the calculations to 
introduce an auxiliary variable @ defined by 


A = A Gen 8, Y = Bsin 6. (61) 


Substitution of these values into Eqs. (53) shows that the coordinates (&, 7) are given by 


£ =a cos 8, 7 = b sin @. (62) 


The inclination ¢ of the normal at (X, Y) to I and the inclination 6, of the normal at 


(£, n) to y are given by 





tan ¢g = 3 tan 6, tan 6, = > tan @. 63) 

The radius of curvature p(¢v) of the curve T is 
po = A*B*/(A° cos ¢ + B sin’ ¢)° 61) 

The value of p; is p + a, and here 
a, = —A(A —a)/(A’ cos’ ¢ + B’ sin’ ¢)'” 65) 
In the inner region, WM, = M, = Mo, and Mo, is zero on y so that 
Q ~9Moa 66) 
on 
1.0 
x/B and y/B 

Fic. 7. Variation of thickness across major and minor semi-axes for ellipse under const ssure 


A = 1.5B). 
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In the annular region, Q. = 0 and Q, is given by Eq. (25). The value of Q, is therefore 


dM oc _ . 
=~, £08 (On — #); (67) 


0), _ Q; cos (6, — ¢) 
evaluated on y. Equating the two values (66) and (67) of Q, leads after some reduction to 


1/2 


(68) 








D(e) ab _—_(A*b® cos’ ¢ + B’a’ sin’ ¢) 
=P a’ + b’ (A’ cos’ g + B’ sin’ ¢)'”*(Ab cos’ ¢ + Basin’ ¢) 


5PPs — * a on 
p 
This equation determines the function D(¢) and substitution in Eq. (60) gives the value 
of Moc in terms of the coordinates (a, ¢). 

Figure 7 shows the variation of the thickness A along the major and minor axes of the 
ellipse for the particular case A = 1.5B. 

The volume of the face sheets in the minimum volume design can be obtained without 
determining the design explicitly. The plate is at collapse under a mode w such that the 
rate of dissipation D per unit area of the plate has the value D = ko,Hh. During collapse 
in the mode w, the pressure p is the only external force which does work and it follows 


that 


7 | pwdA = / DdA = 3koHV. (69) 


where the integrals are taken over the area of the plate and where V is the volume of 
the face sheets. For a circular plate of radius B under constant pressure, the volume 
of the face sheets is given by [4] 
4 
. _ ™pB - 
V = 0.117 +} - (70) 
ool 
Equation (69) was used to compare the average thicknesses of the face sheets of an 
elliptical plate ‘of semi-axes A and B and a circular plate of radius B, both plates 
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Fic. 8 
and B and circular plate of radius B. 
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being subjected to the same constant pressure p. The ratio R of the average thicknesses 
is shown in Fig. 8 for the range 0 < B/A < 1. The volume of the face sheets for the 
elliptical plate is given by 


I 
0.117 2 


7 TAB’. (71) 
To 


II 
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ON THE STEADY-STATE THERMOELASTIC PROBLEM FOR THE 
HALF-SPACE AND THE THICK PLATE* 


BY 


I. N. SNEDDON AND F. J. LOCKETT 
The Department of Mathematics, The University of Glasgow 


1. Introduction. In a recent paper, Sternberg and McDowell [1] discussed the 
problem of determining the steady-state thermal stresses and displacements in a semi- 
infinite elastic medium bounded by a plane. In particular they succeeded in proving 
that the stress field induced by an arbitrary distribution of surface temperatures is 
plane and parallel to the boundary. The problem was treated by the method of the 
Green’s function. 

This paper deals with the determination of the steady-state thermal stresses in 
both a semi-infinite elastic medium and a thick elastic plate. As in the analysis of Stern- 
berg and McDowell, the problem is treated as one in the classical theory of elasticity. 
The method of solution employed is that of double Fourier transforms. A general solution, 
corresponding to an arbitrary temperature field, is obtained in the form of double 
integrals and it is confirmed that the stress field is plane and parallel to the boundary 
of the medium. The particular solution corresponding to axially symmetrical temper- 
ature fields is deduced anda solution found of the problem in which the surface 
temperature is uniform over a circular region of exposure and is zero outside. In this 
special case an expression is given for the difference, | ¢, — a |, of the principal stresses 
in terms of tabulated integrals, and the three-dimensional analogue of the isochromatic 
lines constructed (see Fig. 2). 

The corresponding analysis of the steady-state displacements and stresses produced 
by arbitrary distributions of temperature on the surfaces of a thick plate is then given. 
It is shown that Sternberg and McDowell’s result (that the stress field due to arbitrary 
surface temperatures is plane and parallel to the boundary) holds for a thick plate as 
well as for a semi-infinite solid. This result has been derived by McDowell [2], using 
a method similar to that of the paper [1], and by Muki [3], using a method combining 
the theory of Fourier series and that of Hankel transforms of integral order. The method 
developed here is just as simple as these for deriving the general result, and appears 
to be much more suitable for the discussion of special problems. To illustrate the use 
of the method, the isochromatic surfaces within a thick plate are constructed for a 
special distribution of surface temperatures. (see Fig. 3). 

2. The basic thermoelastic equations of equilibrium. We consider an elastic body 
whose boundaries are parallel to the plane z = 0, and in which is established a temper- 
ature field T + @(x, y, 2), where 7’ is the temperature of the solid in a state of zero stress 
and strain. We assume that there are no body forces within the solid and that its surfaces 
are free from traction. If we take a typical length 7 as our unit of length, the reference 
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temperature 7 as the unit of temperature, and the rigidity modulus yu as unit of stress, 
we find [4] that the equations of thermoelastic equilibrium take the dimensionless forms 

V*u + (6° — 1) grad A = b grad 8, (1) 
where u denotes the displacement vector, A = div u is the dilatation, 6? = 2(1 — »), 
(1 — 2v) (» being Poisson’s ratio), and, in terms of the coefficient of linear expansion 
of the solid, 

b = 21 + v)al’/(1 — 2y). 
The variation of @ throughout the solid is determined by Laplace’s equation 
V’'e = 0, (2) 
(in the absence of heat sources), and, in this system of units, the relation between the 
stress tensor 7,;; and the displacement vector u = (u; , U2 , U3) is given by the Duhamel- 
Neumann equation 
ri; = [(8? — 2) A — 18] 8; + (ui; + U;,:). 


Il 


In this last equation u,;,; denotes du;/dx; with the convention that (2, , x2 , 73) 
(x, Y, 2). 

To solve the equations of thermoelasticity in this form we introduce the double 
Fourier transform 
] 


{*(é, 7, z) = 
_ 2r 


/ [ f(x, y, z) exp [i(éx + ny)] dx dy, (3) 


of each physical quantity f(z, y, z) occurring in the problem. The vector equation (1) 
is then easily seen to be equivalent to the set of three simultaneous ordinary differential 


equations 
(D? — Be? — n?)u* — (8? — lém* — (6? — 1)tt Dw* = —ibto*, 
— (8° — l)énu* + (D? — & — B’n’)v* — (8? — 1)tn Dw* = —ihné*, 
— (8° — 1)it Du* — (8? — lin Do* + (8° D®? — & — 9)u* = 5 DO*, 


and the transform of Eq. (2) is 
(D? — ¢’)o* = 0, 
where we have written ¢ = | (¢? + n’)'” | and used D to denote the differential operator 
d/dz. It can be shown that the solution of this set of simultaneous ordinary differential 

equations may be put in the form 

A* = Ee7** + E’e™’, 
(4) 
u* = (A, + Pézle* + (Al + Ptzje™, 
(A, + Pnze** + (Ai + P'nede™, 


= (A, — iPtze"* + (Af — iP’ tz)e™, 


. 
I 


(5) 


where the constants A, , A, , A; , Af, Aj, Aj, H#, E#’ are arbitrary and P, P’ are given 
in terms of them by the equations 
P = [ibE — (6 — 1)(€A, + 9A2 — tf As) "(8 +1)"; (6) 


P’ = [ibE’ — (8? — 1)(EAl + nAt — if ADE"? + 1)”. (7) 
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The expressions given by Eqs. (4) and (5) satisfy the basic equations of thermo- 
elastic equilibrium in the transformed form. To obtain the solution of the original 
equations we make use of the Fourier inversion theorem. 


fie, y,2) =x ff 1, 1,2) exp [-itee + mp] a da, 


to obtain the expression for a physical quantity in terms of its Fourier transform. The 
solutions (4) and (5) are, of course, suitable only for the discussion of problems in which 
the solid region under consideration is bounded by planes normal to the z-axis. 

3. Solution for the half-space. If we assume that the components of the displacement 
vector, and the temperature each tend to zero as z — ©, the solution for the half-space 


z> 0is 
u* = (A, + &Pz)e™™, v* = (A, + nPze™*, w* = (A; — itPz)e**, (8) 
where A, , Az , A; denote arbitrary constants, 

o* = Keo", (9) 
and P is given by Eq. (6). The Fourier transforms of the stress components 7,, , Ty; 5 %» 
are given then by the equations 
Th, = (—$A, — &P2 + &P — itAs — Pele, 

r*, = (—fA, — ntPz + oP — inAs — nfP2e™, 
* = [—i(8° — 2)(EA, + Pz + nA. + 0°Pz) + B'(—fAs + it?Pz — itP) — DE]e™**. 
If, then, we impose the boundary conditions 
0=G(2,¥), Ter = Te = 0, = O, 
on the plane z = 0, we find that 
E = ¢*(, n), (10) 
where ¢*(£, ») denotes the double Fourier transform of ¢(z, y), and that 
(A, , Az, As) = ibp*[2(8* — 1)g°) "Eg, 2, 28). 
For these values of the arbitrary constants A, , Az , As; it is readily seen, from Eq. (6), 
that P = 0 so that the solution (8) can be written in the form 
u* = ibpte **[2(8? — 1°)", 0, 2f). (11) 
For this solution it is easily shown that 
, = Th, = oF = 0, 
for all positive values of z, showing that the stress field is plane and parallel to the 


boundary in agreement with the result of Sternberg and McDowell. 
If we invert Eq. (11) by the Fourier inversion theorem, we find that the displacement 
vector is given by the equation 
“ d 


b SS ee 
e=pgiy |. | in - exp [-iee + w) - IES 
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In the case in which the prescribed surface temperature (x, y) is axially symmetrical, 
so that ¢(z, y) = ¢(p), where p = (x” + y’)'”, we find that 


o*(, 0) = 2x) [aa | 06(0) exp [itp cos (@ — x)] do, 


where we have written § = ¢ cos x, 7 = {sin x, = pcosd, y = p sin &. We there- 
fore find that, in this case, @*(¢, 7) = *(¢), where ¢*(¢) denotes the zero-order Hankel 
transform of the function ¢(p) defined by the equation 


¢*(f) = | pd o( pf )o(p) dp. 


Denoting the component of the vector u in the p-direction by u, we find that 


Up = uCOSB a v sin 3 


= b[4r(6’ — 1] | o*(f)e** de | i cos (8 — x) exp [—7fp cos (8 — x)] do (13) 


0 
= 4b(8° — 1)" J o*(H)Is(ogde** de. 


-2(4*-1) 
bd “ 

A 

1.0 


O.8- 
O6- 
Z=1 


o.4P 


02 


T 








i i i > 
fs) 2 + 
Fic. 1. Variation of the normal component of the displacement vector in planes parallel to the bound- 
ary of a semi-infinite solid. 
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Also the component in the z-direction is given by 


oe «ie = 7~ [ b*(OJ (pte ® dé. (14) 


J0 


As a special case of the use of these formulae, we consider the problem which is 
considered in some detail by Sternberg and McDowell, namely that in which 


{o> » Sf eo <= i: 
. 0, p>, 


¢(p) = 


in which case ¢*(¢) = J.,(¢)/¢ and, in the notation of [5] we have 
u = 4b¢(8’ — 1) "J(1, 1; —1), w = —43h¢(8’ — 1)’ J(1, 0; —1), (15) 


where 


Je, = [ JAI poe HP dk. (16) 


The integrals J(1, 1; — 1) and J(1, 0; — 1) have been tabulated for ranges of values 
of p and z in Tables 9, 10 and 11 on pp. 545-546 of [5], so that it is an easy matter to 
calculate the components of the displacement vector at any point. For example, Fig. 1 
shows the variation of the normal component, w, of the displacement in planes parallel 
to the boundary. 

Using the formulae developed in [5], it is a simple matter to show that the equations 
(15) are in agreement with the expressions given by Sternberg and McDowell. 

It is of more direct interest to calculate the difference of the principal stresses. In 


our system of units, we have 
ou u 
0, — 059 = 2( Ste ——i. 
Op p 


aX 


o, — oy = b(8" — 1) | SO*()[Jo(ot) — 2i(ek) (os )le™ ae, (17) 
0 


so that, in the general case 


OQ 2 

















Fic. 2. Section by the plane 8 = constant of the surface | ¢, — oy | = constant, in a semi-infinite solid. 
The numbers refer to values of (6? — 1)(¢9 — o,)/bdo. 
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and, in the case in which ¢*(¢) = ¢J,(f), 
go, — ps = bpo(8” — 1)"[J(1, 0; 0) — 2U(1, 1; —1)/p). (18) 


The values of this stress difference at a grid of points in the zp-plane can be calculated 
easily from Tables 6 and 11 of [5], and it is then a simple matter to draw out the lines 
joining points with the same value of this stress difference. The resulting set of curves 
will be curves obtained by cutting surfaces of equal maximum shearing stress by a 
plane # = constant, and will correspond to the isochromatic lines of two-dimensional 
elasticity. These contours are shown in Fig. 2 for the simple problem we have considered 
here. 

4. Solution for a thick plate. The general solutions (5) can be written in a form 
which is more suitable to plate problems, in which the elastic solid is bounded by the 
planes z = + d. The displacements due to the temperature distribution whose trans- 
form is 


6* = E cosh (¢z) + E’ sinh (¢2) (19) 


are given by the equations 
u* = (A, + P’tz) cosh (fz) + (A! + Petz) sinh (¢2), 
v* = (A, + P’nz) cosh (fz) + (Ai + Pnz) sinh (¢2z), (20) 


* = (A; + 2P’tz) sinh (fz) + (Aj + tPfz) cosh (fz), 


=> 
£ 
| 


where 
P’ = (8° — 1)(EAl + nAi + t¢AD/(B’ + De& — tbE’/(8’ + 1)e, (21) 
P = (8° — 1)(€A, + nA, + ifA;)/(8’ + IE — ibE/(B’ + 1)f. 
The components 7,, , T,: , 7, of the stress tensor associated with this displacement field 
are given by the equations 
re = (fA, — t£A; + PE) sinh (fz) + 2Péfz cosh (fz) + ({Al — 1A + P’E) cosh (fz) 
+ 2P’ttz sinh (fz), 


7T*, = ({A, — tnA; + Pn) sinh (fz) + 2Pntz cosh ({z) 
+ (¢{Aj — inAj + P’n) cosh (tz) + 2P’nfz sinh (fz). 
o* = —1(8’ — 2)[(AiE + Agn) cosh (fz) + P¢’z sinh (¢z)] — bE cosh (tz) 


+ 8°A,¢ cosh (¢2) + BitP cosh (tz) + #it’2P sinh (fz) 

— i(8° — 2)[(AfE + An) sinh (tz) + P’t* cosh (¢z)] — bE’ sinh (¢z) 

+ 6’ Ast sinh (tz) + B'i¢P’ sinh (tz) + Bit’zP’ cosh (tz). (22) 
Since the right-hand sides of these equations consist of the sum of an odd and an even 
part, the boundary conditions o, = r,, = T,, = 0, on z = + d, lead to three pairs of 
equations of the type a + b = 0, a — b = O, which imply that a = 0, b = 0. In this 
way we get the three equations 
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(¢{A, — i£tA; + PE) sinh td + 2Péttd cosh fd = 0, 
(tA, — inA; + Pn) sinh td + 2Pntd cosh ¢d = 0, 
—i(8° — 2)[(A,E + Agn) cosh td + Pt’dsinh td] — bE cosh td + 8’ A;¢ cosh ¢d 
+ B'itP cosh td + pi¢*dP sinh td = 0, 
from which, with the help of the second of equations (21), we find that 
(A, , Aa, As) = HbES*(6* — 1)"(E, n, —7f), 


and, in turn from these, that P = 0. 
We get a similar set of equations with the A; replaced by A/ and E replaced by E’. 
In this way we obtain the expressions 


u* = Zibgt-°(6” — 1)~'(E cosh gz + E’ sinh £2), 
v* = fibng “(6° — 1)(E cosh gz + E’ sinh {2), (23) 
w* = $b¢~'(8° — 1)"(E sinh gz + E’ cosh 2), 


for the Fourier transforms of the components of the displacement vector. Substituting 
the values for the constants into the equations (22), we find that r*, = r%& = o* = 0, 
—d<2z<d,so that the Sternberg-McDowell result that the stress field induced by an 
arbitrary distribution of surface temperature is plane and parallel to the boundary holds 
for a thick plate as well as for a semi-infinite solid. 

As in the case of the half-space we can easily derive the solution in the case in which 
the temperature field is axially symmetrical. To illustrate the procedure, we shall con- 
sider the simple situation in which 


6 = ¢(p), on z= +d, 6=0 on z= —d. (24) 
For this distribution of surface temperature, we find that 
a= | s9*(¢) sinh [¢(¢ + d)] cosech (25d) Ju(te) dt, (25) 
v, = 308° — 1 [ 6") sinh [5 + @)] cosech (2rd)Ju(so) dz, (26) 
w = 4066 — 1)" [” 90) cosh [ge + a)] cosech (2rd Jo(Z0) dr. (27) 


Our unit of length (J) is, as yet, unspecified. If we now take 1 = d, so that all lengths 
are measured as ratios of half the thickness of the plate, we find that Eqs. (25) to (27) 
assume the simpler forms 


a= | s9*(¢) sinh [3 + 1)] eosech (2%) Ju(to) de, (25a) 
u, = 46(8° — 1)" [ ow sinh [¢(z + 1)] cosech (2¢)J,(¢p) dé, (25b) 


w = wer - 1" [ " $*(0) cosh [¢(e + 1)] cosech (2¢)Jo(tp) dt. (27a) 
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For any given distribution of temperature on the upper surface of the plate, we can 
calculate the zero-order Hankel transform ¢*(¢) of the temperature ¢(p), and, inserting 
this transform in Eqs. (25a) to (27a), calculate the temperature field and the displace- 
ment field within the plate. In the general case, the evaluation of these integrals would 
be a pretty complicated process, because of the occurrence of the factor cosech (2¢) 
in the integrand. By suitably choosing the function ¢(p) we can, however, obtain 
integrals which can be easily evaluated, and obtain the solution of a representative 
problem. 

For example, if we take the distribution of temperature on the upper surface of the 


plate to be given by the function 


— 4)°6,/8k]{(k — 2)[p? + (k — 2)°)°? — (k + 2)[o + (k + 2)°)°"}, 
(k > 2), (28) 


$(p) = [(k? 


then the temperature at the point z = d, p = 0, or, in our system of units, z = 1, p = 0, 


is 6, , and the zero-order Hankel transform is 


*(p) = [3(k? — 4)76,/k]e~“ sinh (29), (k > 2). (29) 
If, now, we substitute the expression (29) into Eq. (25a) and evaluate the simple integrals 
so obtained, we find that the temperature field is given by the equation 


6 = [(k’ — 4)°0,/8k]{(k —-1—2l[_po + (k -—1-—2)°)7” 


2 ) snare 30) 
—Kk+1+2[P +(kh+142)7°"}. | 











| 





Fic. 3. Section by a plane 3 = constant of the surfaces | ¢, — a» | = constant in a thick plate. The num- 


bers refer to the values of 8k(8? — 1)(03 — o,)/bO0o(k* — 4)*. 


—— 
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Similarly, if we make the same substitution in Eq. (26a), we find that the radial com- 
ponent of the displacement vector is determined by the equation 
up = (O/p{(kK+1+2[o +kK+14+2'7)"” (31) 
—&-—1-a[e + &-—1-8'T"}, 
with 
0 = (k* — 4)°b0,/[16k(8" — 1)]. 
From Eq. (31) we in turn deduce that the difference of the principal stresses is 
o, — og = 20{(k —1 —2Dlp? + (kK —1 -— 2)” 
—&+1+4+a)[0° +&+1+4+2))™” 


+ Ak -—1—2)p [e+ &-—1-2)'7" ” 
—Ak+1+z)e"le +(k+1 +27}. 
In a similar way, we find from Eq. (27a) that 
w = Off? + (kK -—1 —27)' + [or +R +1427}. (33) 
Figure 3 shows the surface distribution of temperature in the case k = 3, together 


with the sections (by a plane ? = constant) of the isochromatic surfaces | ¢, — a) | = 
constant. If we had chosen a higher value of the parameter k we should have obtained a 
surface temperature distribution which was less highly concentrated in the neighbourhood 
of the point p = 0; on the other hand, if we had chosen a smaller value of k, such as 
k = 2.1, the curve would have been concentrated into a small band surrounding the 


origin. 
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BOOK REVIEWS 


(Continued from p. 130) 


and other representations of the general solution of the displacement differential equations of equi- 
librium (e.g., Boussinesq-Papkovich and Galerkin), as well as a clear exposition of the basic concepts 
of function space methods in elasticity. Chapter 7 deals with variational methods, such as the Rayleigh- 
Ritz procedure, the method of Lagrangian multiplier and application of the function-space methods. 

The thermodynamic approach adopted in Chapter 8 (entitled ‘“Thermoelasticity’’), essentially 
that appropriate to a simple homogeneous substance, falls short of recent developments on the subject 
[see, e.g., M. A. Biot, J. Appl. Phys., 27, 240-253 (1956) and M. Lessen, Quart. Appl. Math. 15, 105-108 
(1957)]. Furthermore, the treatment of the subject as given in Chapter 8 is quite elementary and may 
be found elsewhere [e.g., at various places in Zemansky’s Heat and Thermodynamics, McGraw-Hill, 
1957, and Searle’s Experimental Elasticity, Cambridge, 1908, pp. 20-29]; the statement on p. 159 that 
Gibbs entropy relation “. .. holds for all processes, reversible or not’’ is misleading, as the proof of its 
general validity for irreversible processes has not been established. The remaining two chapter-headings 
are Time-Dependent Problems and Nonlinear Elasticity. Chapter 9 is a cursory treatment of vibrations 
of and wave propagations in isotropic elastic bodies. Finally, Chapter 10, following a short account of 
nonlinear elasticity (about 7 pages), deals with the problem of elastic stability. 

On the whole the book is readable and provides a useful summary of some topics in elasticity. 
Despite its title, however, this book is not on a par with Theoretical Elasticity by Green and Zerna, or 


Sokolnikoff’s Mathematical Theory of Elasticity. 
P. M. NaGuor 


Celestial mechanics. By E. Finlay-Freundlich. Pergamon Press, New York, 1958. viii + 

150 pp. $7.50. 

Celestial mechanics is a very old branch of applied mathematics and one would suppose that very 
little new can be said about it. This is true as far as fundamentals are concerned. But interest in the 
subject has been renewed in recent years partly because of the “space age’’ concern with rockets and 
artificial satellites and even more because of the astronomer’s endeavor to learn more of the nature of 
binary stellar systems. Professor Finlay-Freundlich, who has just retired from his professorship in St. 
Andrews University in Scotland, has here provided a new and attractive brief summary of the basic 
principles of celestial mechanics with some emphasis on current applications. 

The book contains a very clear statement of the n-body problem and the reason why it is unsolvable 
in the general sense (Bruns-Poincaré theorem). This is followed by chapters on the application of the 
Hamilton-Jacobi theory (transformation theory of mechancis) to the three-body problem, and the 
theory of perturbations. The two-body problem for extended deformable bodies constitutes the really 
modern material in the book, and there is a good discussion of close binary systems with tidal deforma- 
tions. Pedagogically speaking the weakest part of the book is the chapter on the application of relati- 
vistic mechanics to the two body problem. Here the relativistic four-dimensional line element is intro- 
duced without explanation and the ensuing derivation of the Einstein formula for the advancement of 
the perihelion of Mercury assumes the form merely of a mathematical exercise. 

On the whole the author is to be congratulated on the clarity of his treatment and the large amount 


of fundamental material he has compressed into so small a space. 
R. B. Linpsay 


Testing statistical hypotheses. By E. L. Lehmann. John Wiley & Sons, Inc., New York, 
and Chapman & Hall, Ltd., London, 1959. xiii + 369 pp. $11.00. 
This volume is based, in part, on a well-known set of lecture notes of the author. These notes also 


included a discussion of the theory of estimation. The book is limited to a discussion of testing problems. 
This may be due to limitations of space and perhaps also a belief that the relevant theory in this domain 


is in a more coherent form. 
(Continued on p. 182) 
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DUALITY IN QUADRATIC PROGRAMMING*}+ 


BY 


W. 8S. DORN** 
New York University 


Abstract. A proof, based on the duality theorem of linear programming, is given 
for a duality theorem for a class of quadratic programs. An illustrative application is 
made in the theory of elastic structures. 

1. Introduction. Recent interest in quadratic programming has resulted in a series 
of computational methods for this type of problem. Some of these are described in 

, 2, 3, 4, 5]. Little emphasis, however, has been placed thus far on the concept of duality 
in quadratic programs. This concept, which has proved so valuable in linear programs, 
is investigated briefly in what follows. 

In Sec. 5 a dual problem for a class of quadratic programming problems is formulated 
and the equality of the two objective forms is verified. Dennis [6] has indicated previously 
that such a duality existed based on the Kuhn-Tucker “equivalence theorem” [7]. The 
proof given here rests on the duality theorem for linear programs. 

2. Notation. In what follows, matrix notation will be employed. Lower case letters, 
x, y, «++ will denote column vectors and capital letters A, C, --- will represent matrices. 
Prime denotes transpose so that x’, y’, --+ are row vectors. The product z’y is the inner 
product of the two vectors x and y. 

A vector inequality will apply to each component of the vector, i.e., x > 0 indicates 
that each component of x is non-negative. 

3. Duality in linear programming. The linear programming problem may be posed 
as follows. To minimize the linear form p’x over all n-dimensional vectors x satisfying 
the constraints 

Az > b, 


xz > 0, 
where p is ann X 1 vector, bis an m X 1 vector and A isan m X n matrix. 
The dual problem to the above is to maximize b’v over all m-dimensional vectors v 
satisfying 
Av <pD, 
v> 0. 


The duality theorem [8, 9] states that if a solution to either problem exists and is 
finite, then a solution to the other problem also exists and indeed 


Minimum p’x = Maximum bb’. . (3.1) 
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4. Aclass of quadratic programs. A class of programs which has received considerable 


attention [3, 5, 6] is 


Minimize: f(x) = 3a’Cx + p’x (4.1) 

subject to 
Az > b, (4.2) 
zs 2 0, (4.3) 


where C is a symmetric, positive semi-definite, m X n matrix and p, b, A and z are as in 
Sec. 3 above. This problem will be referred to as Problem I. 

The symmetry restriction on C results in no loss of generality, while the positive 
semi-definiteness requirement assures that /(x) is convex and that a local minimum 


is also a global one [3, 5). 
In order to prove a duality theorem for this class of programs the following lemma 


is required, 
Lemma. If C is a symmetric, positive semi-definite matrix then for any vectors x and y 


y’Cy — x'Cx > 2x’Cly — 2). 
Proof. From positive semi-definiteness, for any x and y 
(y — x)'C(y — z) => 0, 
y’Cy > 22'Cy — x’Cr. 
Subtracting x’Cz from both sides 


y'Cy — 2’Cxz > 2x'Cly — 2). 


5. A duality theorem for quadratic programs. A dual problem to Problem I is 
Maximize: g(u,v) = —}u’Cu + b’v (5.1) 

subject to 
Av — Cu <p, (5.2) 
> (5.3) 


= =, 


where wu is ann X 1 vector and v isan m X 1 vector. This problem will be referred to as 


Problem IT. 


Theorem (Dual). (i) If x = x, is a solution to Problem I then a solution (u, v) = 
(2%) , Uo) exists to Problem II. (ii) Conversely, if a solution (u, v) = (uo , ve) to Problem 
II exists then a solution which satisfies Cz = Cu, to Problem I also exists. In either case 

Max g(u,v) = Min (f(z). (5.4) 

Proof. (A) Suppose first that x = 2, is the minimizing solution of Problem I. Con- 


sider the following linear programming problem 
Minimize: F(x) = —}a(Caxo + x6(Cx + p’r (5.5) 


subject to 
Ax 


IV 
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Denote this as Problem I’. Notice that the constraint sets for Problem I and Problem I’ 


are identical. 
Now suppose there exists an 2* satisfying the constraints and such that 


F(a*) < F(ao), (5.8) 


(xoC + p’)(a* — x) < 0. 
It is easily verified that 
XZ, = Xo + k(x — Qo), e< es <3 
also satisfies the constraints. Now 
f(a;) — f(x) = kl(xeC + p’)(x* — 20) + 3k(a* — x)’C(x* — X)]. 
Choose k to be 


(oC + p’)(* — 20) 
=s 1(x* — ao) C(x* — 2) 


[t follows that the term in square brackets is negative so 
f(z.) — f(t) < 0. 


But f(xo) < f(x:) since x» is the minimizing solution of Problem I so the inequality 


(5.8) cannot hold for any 2, i.e., for all z 
F(x) > F(x). 


Thus x») minimizes F(x) and is the optimal solution to Problem I’. 
The dual problem to Problem I’ is (Sec. 3) 


Maximize: G(v) = —}ajCrxo + b’v (5.9) 

subject to 
A’v < Crp + p, (5.10) 
v> 0. (5.11) 


» 


Denote this as Problem II’. By the duality theorem for linear programs, (3.1), 
Max G(v) = Min F(x) = F(x). 
If v = is a maximizing solution of Problem II’ then the last equation becomes 
bo = x6CXo + p'Xo - (5.12) 

Consider now admissible solutions (u, v) to Problem II. In particular (x , vo) is 

admissible. Now from (5.1) 
gay » Vo) — g(u,v) = —4axiCxo + b’vo + Zu’Cu — bv 

by the lemma 


G(X , Vo) — g(u,v) > x6C(u — 2%) + bo — bv 
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and from (5.12) 
g(Xo 5 Vo) — g(u,v) > roCu + p’x — bv. (5.13) 
Now from (5.2) and (4.3) 
xiCu > «(Av — p) 
and from (4.2) and (5.3) 
—b’v > —25A"v. 

Substituting these last two inequalities in (5.13) 

g(Xo ,%0) — glu, v) > 2A” — cop + p’X%o — %A’v = O. 
Thus (2% , vo) maximizes Problem II. Finally from (5.1), (5.12) and (4.1) 

G(Xo , Vo.) = —FxsCro + bo = 426Cro + p’x%o = f(Xo) (5.14) 


which verifies the equality of the objective functions (4.1) and (5.1). This completes 
the proof of part (i) of the theorem. 

(B) The converse will be proved by applying the above result to Problem II. Suppose 
a maximizing solution (uo , vo) of Problem II exists. Problem II may be rephrased 


Minimize: —g(u,v) = 3u’Cu — b’v 


subject to 
Cu-—- Av> —p, 


v> 0. 
Now let 
u=fr-—8 
where 
r20, 
s> 0. 
Problem II then becomes 
[ ¢c. —C, Ol fr) r 


Minimize: —G(r,s8,v) = 3(r,s,v)’|—C, C, O} js} + (0,0, —b) 








subject to 
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This is now in the form of Problem I, and by part (i) of the theorem (already proved 
in A above) implies the existence of a solution to a dual problem which is 
( c, -C, 0) fw 


Maximize: —}(w, y, z)’|—-C, C, Ol ly| — p’x 
| 0, 0, Ole 
subject to 
Cx — Cw + Cy + 02 < O (5.15) 
—Cxr + Cw — Cy + 02 < 0 (5.16) 
—Azx+ Ow + Oy + 02 < —b (5.17) 
z>0o0 
Moreover, the maximizing solution is required to satisfy 
w—- Y= UW, (5.18) 
z=v 
Inequalities (5.15) and (5.16) imply that 
Cx = C(w — y) (5.19) 


so the dual problem may be rewritten 
Maximize: —}32’Cx — p’x = —f(x) 
subject to 
Az > b, 
z > @, 
which is exactly the original Problem I. From (5.18) and (5.19) then the optimizing 
solution x to Problem I must satisfy 
Cz = Cm. 
Finally from Eq. (5.14) it follows that 
Min —g(u,v) = Max —f(x), 
which completes the proof of part (ii). 
6. Computation of the dual variables. The proof in the preceding section also pro- 
vides a means for calculating the dual variables (uo , vo) once the primal variables, 2» , 


have been found. 
The vector up is identical with zx, . The vector vp is then a solution of a linear pro- 


gramming problem 


Maximize b’v 


subject to 
A’v < Cx + p, 


v> 0. 
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7. Other classes of problems. The quadratic problem (Problem I) may be formulated 
in various other ways with resulting changes in its dual (Problem IT). Some of these, 
including the original, are tabulated below for reference. 


Primal Problem 


Dual Problem 








Type I 





Min 3 2’Cx + p’x 
Az >b 
z > @ 


Max — 3 u’Cu + b’v 
Atv —Cu<p 
v>0 








Type II 





Min 3 2’Cx + p’x 
Ax >b 


Max — 3 u’Cu + b’v 
Ay — Cu=p 
v>0 





Type III 


Max — 3 u’Cu + b’v 





Min 3 2’Cxz + p’x 








Az = b Av—Cu<sp 
soe 
Type IV 
Min 3 2’Cx + p’x Max — 3 u’Cu + b’v 


Ar =b A’vy — Cu =p 

The Type IV problem may, of course, be treated by standard Lagrange multiplier 
techniques. The dual problem for a problem of this type has been given previously 
[10]. Indeed, v are the multipliers for the original problem. 

Notice that at the optimum, in all types listed above, 
(7.1) 


Uu z. 


8. An application to elasticity. As an application of the duality theorem for quad- 
ratic programs, consider the problem of determining the elastic solution of a plane 
pin-jointed truss consisting of n bars and m joints (n > 2m — 3). The truss is externally 
statically determinate and the applied loads lie in the plane of the truss. 

The problem may be formulated as minimizing the strain energy subject to equi- 
librium constraints (Castigliano’s Second Theorem). 

If S; denotes the force in the jth member and A, , £; , L; are its cross-sectional area, 
elastic modulus and length respectively, then the strain energy U is 


raters & 


i=1 ? 
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and the equilibrium conditions may be written [11] 
2 Gull = ff, («= 1,2, --- ,2m — 3), 
i=1 


where F; is the force component at the ith joint. The a,;; depend on the geometrical 
configuration and are essentially direction cosines of the angles between the bars and 
the coordinate axes. 

This is a problem of Type IV and the dual problem to this minimum problem is 


from Sec. 7 


si — Lb «@, eo 
Maximize: —} >- AE Sit Do Fu; 
gaa eos i=1 


subject to 


<— L; ; 
a S. =f me —— * 
d a; ;U; AE, > ¢) 1,2 n) 


The use of the variables S; in the dual is justified by Eq. (7.1). Making use of Hooke’s 


law which gives the elongation, e; , of the jth bar as 


a Lj 8 
” da * 


é 


the problem becomes 


n 2m—3 


ial A,E; 2 , 
Maximize: —} >>  % e+ > Fu, 


subject to 


d au; —e; = 0 (9 = 1,2, +--+ ,m). 


If, for the moment it is assumed that the u, are the displacement components of the 
ith joint, the objective function for the dual problem is the negative of the total potential 
energy and the constraints become the compatibility equations. The dual problem is, 
therefore, equivalent to the Theorem of Minimum Potential Energy. The equality of 
the objective functions results in a restatement of the Principle of Virtual Work. 
Other applications to electrical networks containing resistors, diodes and voltage 


and current sources may be found in [6]. 
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CANONICAL AND HAMILTONIAN FORMALISM APPLIED 
TO THE STURM-LIOUVILLE EQUATION* 


BY 
M. A. BIOT (Shell Development Company, N. Y.) 
AND 


I. TOLSTOY (Columbia University, Hudson Laboratories, Dobbs Ferry, New York) 


Abstract. The Sturm-Liouville equation is expressed in Hamiltonian form. A 
simple generating function is derived which defines a large class of canonical transforma- 
tions and reduces the Sturm-Liouville equation to the solution of a first order equation 
with a single unknown. The method is developed with particular reference to the wave 
equation. The procedure unifies many apparently diverse treatments and leads to new 
insights and procedures. Some new transformations are obtained, useful in the turning 
point region and for the improvement of accuracy in the region of validity of W.K.B. 
solutions. In addition a new power series expansion near the turning point is obtained. 

1. Standard canonical form of the Sturm-Liouville equation. Let q, p be conjugate 

ariables and f, g arbitrary functions of the independent coordinate z. The Hamiltonian: 


1/1. 2 
n=1|> + ot’ | (1.1) 
defines the canonical equations 
dq ] 
lz 
; f (1.2) 
1 
ale 


which are equivalent to the second order equation 


3 (1 “2) Hes (1.3) 
i.e., the general Sturm-Liouville equation, having a vast number of physical applica- 
tions. A very large fraction of the integrable equations of mathematical physics are of 
this type. Thus, in separable coordinate systems, the solution of the partial differential 
equations of Laplace, of heat flow, of mechanical electro-magnetic and quantum- 
mechanical waves reduces to some form of Eq. (1.3). However, the exact behavior of 
the solutions of this equation is known only for a relatively few simple cases. In problems 
of practical interest it occurs frequently that the functional forms f and g are such that 
no solutions of Eq. (1.3) in terms of known functions is possible. One is thus often re- 
duced to devising approximations and/or numerical schemes of integration. 

Hamilton’s equations (1.2) are a system of first-order equations entirely equivalent 
to the Sturm-Liouville equation (1.3). From a purely numerical standpoint one might 
expect that more advantageous formulations are possible and, since we are dealing with 


*Received August 17, 1959. 
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a set of canonical equations, it seems only natural to apply the theory of canonical 
transformations which provides a natural and flexible means of exploring equivalent 
systems of first order equations. In the following sections we will describe one group of 
transformations, which should not only have certain numerical advantages, but also 
provide new insights and results of a general nature. 

2. Canonical transformation of the standard form. Clearly, considerable advantages 
can accrue from finding an equivalent set of first order equations, one of which is un- 
coupled, the other being soluble by quadrature. These requirements may be satisfied 
if one defines equivalent canonical equations in the new variables corresponding to the 
following class of generating functions: 


F(q, P,z) = 7 ®(P,2), (2.1) 


where @ is an arbitrary function. The corresponding transformation is: 


oF 
= -= 9 9 
p aq 2q® (2.2) 
oF 2 Ob 
—- — — .—— 2 e 
@ ap 7 ap ae 
and 
, Q 
" = — 9 
q ag (2.4) 
oP 
1 _ Qe 
p= ab (2.5) 
oP 


The new Hamiltonian: 


ae ae oF i111. é » O® 
- ad 2 2 > ai. oe = 2 ! 2 _ specif 9 > 
H=>5 E p + 949 | et E p+ 9q | tT? (2.6) 


becomes, after substituting (2.4), (2.5): 


H = QK(P, 2) (2.7) 
with 
5 OP 
zi 1| 46° g | Oz : 
(P,2z) = =| a . (2. 
aisle 2|,db db! | ae ain 
|? oP oP oP 


The canonical equations for Q, P are 


Te ss ll (2.9) 
dz 
dQ will OK ; (2.10) 


dz oP 
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The last is integrable by quadrature: 


"OK 
Q = Q, exp i 9P ae| (2.11) 
3. Selection of transformations. Assume that 
¢(P, 2) = o(2)W(P) (3.1) 
then 
d 
, Log ¢ 
- 1) 4e¥ = 2 dz — ae 
ies 2 dv 7 dy + d | e (3.2) 
lap %ap| ape 
One may write, without loss of generality: 
S(P en 
ae (3.3) 


¥ = RP) 
This facilitates the selection of the Q, P that are most desirable for a given problem or 


class of problems, since now 


q= Q’’RP) g'? A (3.4) 
p = 2Q'?S(P) ¢? a”, (3.5) 
with 
dS - dR ‘ 
A —_ R dP — dP (3.6) 
In this case: 
, l _ 2 _d ™ 
K=5 {1 F gS? + oR | + RS |. Log e}. (3.7) 


Equations (3.4), (3.5) enable one to select explicitly the forms R(P), S(P) in an appro- 
priate fashion, since R, S are arbitrary. This approach gives a free rein to physical 
intuition. Generally speaking, in the selection of transformations one may adopt one 
of two attitudes. 

On one hand, one can look for special properties that would be advantageous for a 
specific problem. Thus, if one is in a position to anticipate good formal approximations 
to g, p the functions R, S, can be chosen so as to produce the corresponding forms in 


Eqs. (3.4), (3.5). This will result in almost linear or almost constant solutions to Eqs. 
2.9), (2.10) and thus greatly facilitate numerical integration. But this approach is 


fruitful in a relatively small number of cases and will not be discussed here. Besides the 
various perturbation schemes, such as that of Miller and Good [1] are well suited to 
this kind of situation. 

On the other hand, one may wish for properties suitable for integration in many 
cases, so as to provide a generally useful method; this is essentially the point of view 
adopted in the following sections, and which we will illustrate in the context of the wave 


equation. 
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Note also that the choice of transformations need not be based upon Eqs. (3.4), 
(3.5), but can stem from formal considerations concerning the desired behavior or form 
of individual terms in K(P, z), Eq. (3.2). This approach is less intuitive than the one we 
have followed. We shall have occasion to remark further upon some of these transforma- 


tions. 
4. The wave equation: a polar type of transformation. If in Eq. (1.3) we take 


f=1, 
(4.1) 
g=B =h(2) -Fk, 
where k is a parameter, we have the time-independent wave equation 


dq 
dz” 


+ Bq = 0. (4.2) 
This equation is typical of equations with oscillatory solutions and thus provides parti- 
cular motivation for the use of transformations of the class discussed in the preceding 
section. For example, one may look for functions R(P), S(P), satisfying the following 
criteria: 

(1) The unknowns P, Q are to be reasonably monotonic while representing at the 
same time functions qg, p that may be highly oscillatory. 

(2) Neither of these functions, nor any of the coefficients of their defining equations 
shall be singular in the region of interest. 

Criterion (1) is readily fulfilled by taking for R, S forms such that their dependence 
upon P is of the same general type as the expected dependence of q on z. This can be 
achieved approximately, in a broad and general fashion, by taking for FP, S trigonometric, 
exponential or Bessel functions of P. 

Criteria (2) require only that R, S, ¢ be functions without singularities and that 
A, ¢ have no zeros. 

A simple and obvious series of choices meet 


¢ = const = 2”. (4.3) 


these requirements. First of all, let 


Note that this corresponds to dropping the last term in A(P, z), Eq. (3.2), thus greatly 
simplifying the formal problem. 
Next, take 


R = sin a 
(4.4) 
S —cos P, 
this gives 
q = (2Q)'” sin P, : 
f . 1.2 
p = ~ 2Q)' - cos P. 
And by virture of Eqs. 2.9), (3.7), P obeys the equation 
dP —e 
+ cos P + 8 sin’ P 0 1.6) 


dz 





tw 
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and 


Q = Q, exp i (3? — 1) sin 2P az|. (4.7) 


QP, and the coefficients of the differential equation (4.6) are well behaved everywhere, 
including the region of the turning point 6° = 0. Boundary conditions at a discontinuity 
imply simply continuity of P. This formulation has therefore definite merits of simplicity 
in appearance and in analytical behavior, as well as generality. In this sense, it has 
obvious advantages over other, better known, equivalent systems of first order equations 
sometimes used in practice, such as the Ricatti equations for the reflection coefficient or 
for the impedance: the former has oscillatory solutions and a coefficient that is singular 
at turning points, and the latter has solutions that are both oscillatory and singular. 

The transformation (4.5) is, of course, a particular case of the more general group for 
which ¢ is not constant. Because of their aspect, it seems natural to call them polar- 
type transformations. The corresponding generating function 


F = —q ¢ cotanP (4.8) 


is the generalization to time-dependent systems of the usual function for oscillating 
systems found in many textbooks. Obviously, different choices of g may result in trans- 
formations that would be better adapted to specific cases. But, of the transformations 
satisfying the second criterion of analyticity, Eqs. (4.5) give one of the more general 
types since, because of its relatively non-committal form, it has a wide range of appli- 
cation. 

5. A polar transformation for slowly variable parameters: Relation to the W. K. B. 
approximation. It may sometimes be expedient to ignore criteria (2) of the preceding 
section and to permit the appearance of singular coefficients, providing that the singu- 
larities appear outside of the range of intended applications of the resulting equations. 
The following purely physical reasoning leads to a useful transformation of that type. 

One may interpret Eq. (4.2) as the equation for a pendulum of variable length, of 
period 27/8. The average energy per cycle is 

E = Bq’. (5.1) 

If 3 varies very slowly one has, by Ehrenfest’s theorem for slowly varying constraints, 
the adiabatic invariant 

TT 

B 


E = const. (5.2) 


And, since 8 varies slowly, for each cycle: 
- 1 , 
q = F xX const. (5.3) 


The constant term arises from the time average over a cycle of an oscillating, almost 
periodic function, i.e., g is of the form 
-1/24/ - 
es 8 f(2), (5.4) 
where f(z) is an almost periodic function of period 27/8. It is therefore only natural to 
choose ¢ proportional to 6 in Eqs. (3.4), (3.5) e.g., 


g = 36. (5.5) 
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Using script capitals for our new canonical variables Q @ to distinguish them from the 
Q, P of Eqs. (4.5), we have 
aa)" .. 
q=\F> sin @, eG 
p (9.0) 


—(208)'” cos @, 


II 


P 


corresponding to the generating function 


F = —}qB6 cotan @. (5.7) 
This can also be obtained by imposing formal requirements on K(P, z) in Eq. (3.2) 
Thus, with f = 1 and Eq. (5.5) one has 
K = +8 ean + Se: Log 8 5.8) 
— ae — — Log B. (5.1 
dy dy dydz = *" : 
de de de? 


Choosing the — sign in front of the bracket leads to Eq. (5.7) if one imposes the condition 
that the term in brackets be equal to —1. The equations for ®, Q are 


de Ss 2 . x O 
Ps _ 5) q, 08 8) sin 20 + 6 = 0, (9.9) 
Q = Q exp | | cos 20 d Log a |. (5.10) 
An instructive change of independent variable is 
= -| B dz, (5.11) 
giving 
dP ] d ‘ ~ 
— == < 7? y ~ 2°? — = (? ¢ 
ae 2 ds 18 8) sin 26 l 0. 5.12) 
In regions of very slowly varying 6 
d er 
j Log 6B <1 (5.13) 
as 
an approximate solution of Eq. (5.12) is 
O(s) Ys +H, (5.14) 


i.e., ®(s) becomes a straight line of slope unity. Condition (5.13) can be interpreted 
physically: it implies that the change in parameters is so small that the backscattering 
or continuous process of partial reflection of a wave train can be neglected [2]. This 
backscatter and multiple scatter is measured by 


1B . 


d i 
oy {[ = 
ds "® B ds 


It will be small if either d8/ds is small or if 8 is large, (small wave length). Elsewhere & 
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: 











—_ 


Fic. 1. Comparison of exact and approximate solutions of Eq. (5.12). 


will generally be a monotonic wave, gently oscillating about the straight line (5.14), 
as shown in Fig. 1. 


Q => Q = const (5.15) 


and 
q = As’ sin (s + ®), (5.16) 


where A is a constant. This is the first-order W.K.B. approximation. 

We thus have a new interpretation of the W.K.B. formulae: the phase term is the 
limiting form for slowly variable parameters of our canonical @ variable, whereas the 
constant of integration A is the limit of (2@)'”? under the same conditions. 

In the vicinity of a turning point 8 = 0 the criterion (5.13) is no longer fulfilled: 
the variation of parameters per unit z wave length becomes very fast and the reflected 
components become large (since one is in the region of glancing incidence and total 
reflection). Indeed, the coefficient: 


d 1 dh 1 (5.17) 


ds CB = 5G, BF 


is divergent. As could have been foreseen: the representation (5.6) is not felicitous for 
rapidly varying parameters since we had deliberately designed it for slowly varying 8. 
Here, the transformation of Sec. 4 is much more appropriate. It is connected to our 


P variable thus: 


I - 
P = tan'' | —tan@ |. (5.18) 
B 
Che behavior of @ and P far from a turning point is similar only for 8 = 1. Otherwise P 
consists of a sequence of steps crossing back and forth through the @® curve (Fig. 2). 








z 
a i_- 


Fic. 2. General behavior of variables P and &. 
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6. A new turning-point expansion. An expansion of P, Eq. (4.6), in power series 
of 8° is easily achieved, leading to a general result valid on both sides of the turning 
point 6° = 0, for any functional dependence of 8° on z. 

For example, we may write 


tanP=a+a6 +a8 +-::-, (6.1) 
where dy is an integration constant. Substitution in Eq. (4.6) yields, writing dh/dz = h’: 
l 
Se = 
h 
I j ‘ 
ay —— (a + a2) 6.2) 
2h 
— = (2a,a, a 5Q;) . 
on: 
oli 


As pointed out below tan P is the impedance. This gives a physical interpretation to 
the expansion (6.1). It is also easily verified that tan P satisfies a Ricatti equation. 
Applying Eqs. (4.5) and (4.7) and neglecting terms of the order of 8° and higher gives: 

l 


2 l 3,0 s - 
q A] ao — —~B +548 | exp] a 5 + (8) |. ().3) 


he ~ ~ 
To extend this to terms of the order of a}8* and higher, take 


er Po 
(8°) = —7 408 — A | B° dz. 6.4) 
Additional terms are easily secured. These expansions are convergent In the vicinity of 
; ; 
8 = 0, but become divergent in the region P > tan” 1/8. We have not made a rigorous 
study of the radius of convergence of the series (6.1). 
For a small region near a simple turning point (h° # 0) we can approximate 6 by 


a linear function of z — 2 , 2) being the turning-point coordinate. Equation (6.5) then 
gives the wave function as a constant plus a term linear in (z — z,). This same result is 


obtained by taking Langer’s approximation [3], which is of the type: 
y 2 — 20) [BJ i3{(2 — 20) } + CI rah (2 — 20)" $);, 6.5) 


expanding it in powers of (2 — z)) and keeping the resulting constant and linear term. 
Our results (6.1)-(6.3) are more general since they are the correct expansions for any 
functional dependence of 8 upon z, whereas Langer’s formulae are based upon the approxi- 
mation that 8° is linear in a small region surrounding the turning point. 

7. General remarks. ‘The discussion of Secs. 4, 5, 6 would not be complete without 
some brief remarks concerning the physical meaning and previous derivations and uses 
of the P and @ variables in the context of the wave equation. 


The meaning of P is obvious from the defining equations (4.5): 


f 
i ¢ tan P 7.1) 
q 
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is an impedance. The use of P is therefore equivalent to replacing impedances ¢ by 
tan ' ¢ thus eliminating the singularities present in impedance solutions. 

As for ®, a simple physical interpretation in the case of total reflection can be obtained 
as follows. Consider the well-known Ricatti equation for the reflection coefficient ® [14] 
d 


a Log B(1 — &’) + 278R = 0 (7.2) 


dR ©. 1 
dz 2 


the substitution 


R=e" (7.3) 


shows that ¢ obeys Eq. (5.9). In other words, when ¢ is real (total reflection), it is iden- 
tical to ®, which is thus the half-change in phase of a plane, totally reflected wave. 

The transformation (7.3) has been used by a number of authors, notably Eckersley [4], 
Hartree [5], Bremmer [6], Brekhovskih [7] and others for approximate descriptions of 
wave propagation in stratified media. Walker and Wax [8] suggested the use of Eqs. 
(7.3) and (5.9) for the numerical integration of transmission line problems. Tolstoy [9] 
has used exact recursive relations for @ in layered media when Log 8 is a sum of step 
functions, relations that are very similar to those used by Biot [10] in studying the 
torsional oscillations of loaded rods. 

A transformation which is not canonical but very similar to (5.6) 


q = psin 6, 7 4) 


q = p8 cos 8, 


was first introduced by Priifer [11], and used again later by Barrett [12] in investigating 
asymptotie properties of Sturm-Liouville operations. Atkinson [13] pointed out its 
possible applications to numerical problems. Its usefulness in that sense would be about 
equivalent to (5.6). 

8. Conclusions. We have presented a class of canonical transformations which 
reduces the Sturm-Liouville equation to two first order equations having a number of 
interesting properties: an equation in P which is independent of Q, and an equation in 
Q which is soluble by quadrature. The general form of these transformations is quite 
flexible and involves arbitrary functions of P, z. In particular it is easy to generate trans- 
formations of the dependent variable into any desired functional forms, e.g., forms 
having properties considered suitable for a given problem or class of problems. This 
leads to new insights and forms useful for general numerical integration schemes. 

We have contented ourselves here with illustrating the potential of this approach 
by simple examples in the context of the wave-equation. In the process we have been 
led to several interesting results: first, a new transformation (Sec. 4) suitable for numerical 
integration for all values of the parameters, valid on both sides of a turning point. A 
second result is a transformation mostly useful for slowly varying parameters; when this 
variation is slow enough it leads to the W.K.B. results, thus giving a new interpretation 
of these formulae in terms of canonical variables. Thirdly, a power series expansion in 
terms of the wave number 8 leads to a novel and general formula for the wave function 
near, and to both sides of, a turning point, valid for any functional dependance of the 


coordinate. 








~ 
“I 
bo 


Hm CO bo 
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VIBRATION OF QUARTZ CRYSTAL PLATES 


BY 
R. P. JERRARD* 
University of Illinois 


1. Introduction. Among the most difficult vibration problems are those arising 
in the analysis of quartz plate vibration. In the first place the quartz crystal is highly 
anisotropic. Secondly, the most pressing problems involve plates vibrating in complicated 
modes where simple plate theory does not apply. Finally, many of these plates are of 
non-uniform thickness, being bevelled at the edges or ground into convex form. At 
present, the best guides to an understanding are the work of Mindlin and Newman 
[1], [2], Sykes [3], and Ekstein [4]. 

In particular, Mindlin has developed a theory of elasticity for plates based on ex- 
pansion of the displacement functions in powers of the distance from the center plane 
of the plate. He has successfully used this theory in describing the interaction of the 
fundamental thickness shear mode with high harmonics of the lengthwise flexure modes. 


Y 











h (x) 
x= A/2 





ps 


Fia. 1. Coordinate axes and shape of crystal plate. 


teferring to Fig. 1, the crystal is said to vibrate in fundamental thickness shear 
when the top and bottom surfaces are moving in opposition to one another, and their 
motion is predominantly in the x-direction. A lengthwise flexure mode is one in which 
the middle surface of the crystal has a y-displacement.like sin mx. Interaction between 
the two modes takes place only when their frequencies are about equal, and since the 
crystal is quite thin, this will occur only for high harmonics of the flexure mode. 

Very little analysis of plates with non-uniform thickness has been made, though 
many such plates are manufactured. It is found empirically that with a plate of convex 
shape the desired thickness shear modes are purer, that is, freer of distortion which is 
apparently due to coupling with other modes, particularly the high order lengthwise 

*Received April 17, 1959; revised manuscript received September 11, 1959. This work was carried 
out while the author was a staff member at Bell Telephone Laboratories, Whippany, N. J. 








174 R. P. JERRARD [Vol. XVIII, No.2 


flexure modes. It is also found that with a convex, or “contoured” plate the surface 
displacement near the edges is very small compared to the displacement at the center. 

In this paper, an analysis of such contoured plates is made. It is found that with 
one further approximation, Mindlin’s plate theory can be applied to a certain class of 
contoured plates. These solutions have been obtained, and the frequencies computed 
for certain plates. The analysis is given in Sec. 2, and the numerical results in Sec. 3. 

2. Analysis of contoured plate vibration. We take as a starting point the plate 
theory developed by Mindlin [2]. The coordinate system is shown in Fig. 1. The I. R. FE. 
standard notation is used for the stresses, displacements, and the elastic constants of 
quartz. 

There are two essential features of the solution. First, it is assumed that 7. , the 
stress in the y-direction, is zero. This is reasonable for thin plates. Second, the form of 
the displacement components is chosen to be 

u = yW(x,2, U), 

v = nr, 2, 0b), 

w= y¢(r,2z,), 
where ¥, n, and ¢ are functions to be determined. In the expansion of the displacements 
in powers of y we use at most only the linear term. We later specialize to the case where 
w vanishes. It is seen that the y dependence of the displacements has been prescribed so 
as to permit a thickness-shear type of motion. 

The general procedure for obtaining equations of motion is as follows. There are six 
elastic stress-strain relations, one of which is altered when T, is set equal to zero. This 
equation is solved for dv/dy, and this quantity is eliminated from the remaining five 
equations. The assumed displacement functions of Eq. (1) are substituted into the five 
equations, which are then integrated with respect to y. This integration introduces 
three stress couples and two stress resultants which are defined as follows: 


M.,M.,M..)=[ (T,,T:,Tsy dy, 


t 
bo 


(QO. .Q.) = | (7 «5 1) GQ. 


The five stress-strain relations now express these five plate stress components in terms 


of the functions Y, n, ¢ Ol Kq. (1) as follows: 


“.=D,—+D,—, 
or 02 
ry Vee 

M,=D,*+D,—+ 
Or O02 


o 
Q, = D,| ¥ ry in|, 
or 


or 
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“iis 
D, = (Ci, — Cie/Cr2)h*/12, 
Dz = [Cis — Cy2Coa/C22)h*/12, 
Ds = [C33 — C23/C22)h*/12, (4) 
D, = [Cu — C3,/Cr2)he’/12, 
D, = C,,h' /12, 
Dg = Coghn’/12. 


The multiplier 7°/12 in D, and D, was introduced by Mindlin as a correction factor. 
The constants C;; are the elastic constants of quartz. 

We must also make use of the three stress equations of motion. By integrating them 
with respect to y, one can write these equations in terms of the five plate stress com- 
ponents given above. In carrying out this integration, and in the previous integration, 
it is necessary to apply the boundary conditions 7, = 7, = 0 on y = + h/2. These 
conditions simply state that the shear stresses vanish on the top and bottom of the 
plate. Actually, 7’, and 7’, are the shear stresses on a surface normal to the y-axis, and the 
contoured surfaces are not normal to the y-axis. The error involved is small if the con- 
touring is limited. 

The equations of motion in terms of plate stress components can be written fairly 


simply when a factor e’** is introduced. They are 
aM oM,, 
~ = Q. + ie - rz + phiw yp 12 _ 0. 
OX O02 
aq a). > = 
aie a wa + phwn = 0, (5) 
Ox Oz 
aM, aM,. - 2 
- *-Q,4+° 4 + pliwe/l2 = 0. 
Oz Or 


We now assume that the motion of the plate is independent of z. This would be the 
situation if, for instance, the plate were infinitely long in the z direction. With this 
assumption, we find that ¢, W/,, , and Q, all vanish, so the third of Eqs. (5) is auto- 
matically satisfied. The only remaining independent variable is x, and we denote differ- 
entiations with respect to « by primes. Equations (5) reduce to the following two equa- 


tions 


(DW)! — DAv + 7’) + ph’o’y/12 = 0, “ 
(6 
[Dy(y + »’)|’ + pha n = 0. 


These are the equations that must be solved, subject to certain boundary conditions 
at the edges of the plate, x = + //2. 

At this point Mindlin and Forray [5] obtained an approximate solution by putting 
D 0, and introducing a correction factor for D, . Their solution approximates the 
frequencies of the “uncoupled” shear and flexure modes. In the present paper, another 
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solution is obtained which retains the coupling between modes. This solution requires 
that the approximation 
h” oe hh’, (7) 
Pr ‘““ . . 9 . 3 ° ‘ fa 

where / is a constant “effective” thickness, be used wherever h* appears in Eqs. (6). 
The effect of this is to alter the stiffness, bringing it closer to the stiffness of a flat plate. 
In other words, this approximation will reduce the calculated effect of contouring. 

We proceed to carry out the differentiations indicated in Eqs. (6), first introducing 
the notation 


h® 
Din 12 = J), 
(8) 
Deh = Ds , 


so that D,» and Dso do not depend upon the thickness. The equations of motion (6) can 
now be written 


+ = 0, 


» hy ,, 12 D0o¥+ 0" , pw 
} as Po soanet a io a 
¥ + h? p D.. h* ve (9) 


Se ; ” 
Vita" +7 etn) + =*=0. 
h Pes 
This pair of equations can be reduced to a single equation by first introducing the approxi- 
mation (7), and second assuming a solution of (9) of the form 
y = (co — 1)’, (10) 


where ¢ is a constant. 
After using (7) and eliminating y by means of (10), Eqs. (9) become 


h’ pu eg 12 De 
ee _ 3° In’ = 0, 
m . * E o— ] DK h 


— eee od (*'Y we 
7 [" “tare” wre 


where the second of Eqs. (9) has been differentiated. 
These two equations can now be made identical by adjusting the constant o. The 
value of ¢ which makes them identical is obtained by equating the third terms of the 


(11) 


two equations to each other 


pw o L? 2. pu (““) 
! = a os : 12) 
“Ssyar * 2a. a —_ 


We note that if o is to be a constant, then the right hand term in (12) must be constant. 
This restricts the form of h, for then we must have 


hh’ , 
( ) constant. (13) 


This is a differential equation in h(x), which can easily be solved to yield 


h(x) = he exp (i. a + kat). (14) 
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When h(x) is any expression of the form (14), the approximate solution can be carried 
through with no further assumptions. We will consider only the family of crystal shapes 
that result when k, = 0, k, < 0. 

Returning to the constant ¢, one can solve the quadratic equation (12) for o to obtain 


1 7 1 h’ , i h’ , 32 aa)" , 
, sito -a(P) + ([1 0 +(e) R \, (15) 


where 
_ Doo 
g —* Dre ? 
A= Doo , (16) 
pw 
h? 
R = —- 
: 12 


Note that the frequency occurs in the constants J, ¢. 
We now have one solution of Eqs. (9), which is given by Eqs. (10) and (15). The 
remaining equation can be taken as the second of (9), which is now written in the form 


h 


nota fet gud (17) 
n Ao 
There are two values of 1/¢ given by Eq. (15), and since (17) is of the second order, 
there will be two solutions of (17) for each value of ¢. The exact form of (17) depends 
upon the form chosen for the thickness h(7). 

The boundary conditions are given by Mindlin [2]; the edges x = + 1/2 are free, 


and the boundary conditions are 
M, = Q. = 0, x= +1/2. (18) 


Solutions of (17) must satisfy these conditions. Using Eq. (3), the conditions can be 


rewritten in the form 


Mia at z= +1/2. (19) 
¥+7'= oJ 
Thus, corresponding to the four solutions of Eq. (17) (two for each a), there are four 
boundary conditions to determine the arbitrary constants. 
We can now indicate the complete solution to the problem. Taking h(x) as given 
by Eq. (14), the differential equation (17) becomes 


n’’ + (kit + ke)n’ + BE 7 = 0. (20) 
\o 


For each ¢ this equation will have two solutions, one even in x, the other odd in x. Since 
the fundamental thickness shear mode of vibration will couple only with the even 
flexure modes, we consider only the solutions to (17) which correspond to even flexural 
modes. For these modes, 7 is an odd function of x, so our solutions will be of the form 


5; = A,f(o; » z) + Aof(a2 » 2), (21) 
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where f(c, x) is the odd solution to Eq. (20). We need apply the boundary conditions 
only at x = 1/2, for since the solution is odd in 2, the conditions will automatically be 
satisfied at the other edge. 

From (21) and (10), the form of y will be 


Y = A,(o, — 1)f'(o, , x) + A2(o2. — I)f’(o2 , 2). (22) 


Substituting these expressions for ¥ and 7 into the boundary conditions (19), one gets 
two linear, homogeneous algebraic equations in the two unknown constants A, and A, . 
To obtain an equation for the frequency, one sets the determinant of the coefficients in 
the two equations to zero. This latter condition, which provides the frequency equation, 
must be satisfied if the boundary condition equations are to have a solution. Proceeding 
in this way, one gets the frequency equation 


oi(a2 — 1)f’’(o2 , 1/2)f’"(o, , U/2) — on(o, — If’’(o, , l/2)f’(o2 , 1/2) = O. (23) 


The only unknown in this equation is the frequency w, which appears in both o;, 
and o, . In general, Eq. (23) is a transcendental equation, and will have an infinite 
number of roots. This is to be expected, since any crystal of fixed dimensions can theo- 
retically vibrate in an infinite number of overtone modes. Thus for any given set of 
dimensions, one can solve Eq. (23) to obtain the natural frequencies of vibration. 

3. Numerical results—conclusions. In order to obtain the frequency spectrum 
for this class of quartz plates, we must discover the distribution of roots of the trans- 
cendental equation (23). This equation involves the function f(¢, x) which is the odd 
solution of the differential equation (20). The solutions of this equation are certain 
confluent hypergeometric functions [6], which are not available in tabular form. Thus 
the task of finding the roots of (23) involves a very lengthy calculation. 

The expression on the left side of (23) depends, aside from elastic constants, upon 
a total of three dimensionless parameters; when two of these are specified, one can 
obtain a sequence of values of the third parameter at which (23) has a root. This must 
be done for a large number of combinations of the first two parameters in order to 
obtain a good picture of the frequency spectrum. An IBM 704 computer was used to 
carry out this calculation. The required functions were expressed as series, and a pro- 
cedure similar to Newton’s method was used to find the roots. 

The left side of (23) can be expressed in terms of the four parameters 


w/a, U/h, k,l’/8, g- (24) 


Here / is the length of the crystal plate in the x-direction, g(defined in Eq. (16)) is equal 
to 0.283 for AT-cut quartz, h is an average thickness as defined below, and w is the 
radian frequency of vibration. Also @ is the frequency in fundamental thickness shear 
vibration of an infinite flat plate of thickness h. We see from Eq. (14) that 

k, 0/8 = log [h(l/2)/ho] = logr, (25) 


where 7 is the ratio of edge thickness to center thickness. 
The effective thickness h will be defined by 


1/2 1/2 
h = |? | h*(x) az | , (26) 


J0 
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This definition is obtained by minimizing the integrated squared difference between 
h*(x) and h*, since we are approximating the former by the latter. For h(x) as given 
by Eq. (14) we find (4, < 0) 


h = ho(—2/2k,2)"*lerf(—k,2/2)'?)'”. (27) 


Experimental data would show whether definition (26) or some other weighted average 
would be better. However, the definition of h is not used in the computations described 
below. 

The calculation of frequencies was carried out for five values of r, namely 1.0, 0.778, 
0.606, 0.448, and 0.368. That is, for each of the values of r, the frequency ratio w @ was 
obtained as a function of l/h. Figure 2 shows the resulting frequency spectrum when 
l/h is in the neighborhood of twenty-eight. 

When log r is zero (r = 1), the curve coincides with Mindlin’s result [2] for the 
flat plate. The curves have the characteristic structure that also appears in the flat 
plate solution. The more nearly horizontal portions correspond to vibration which is 
predominantly thickness shear, while the more vertical portions correspond to flexural 
vibration. This analysis does not show any marked change in the interaction between 
these modes as the contouring is increased; the main effect is to bring the frequency 
ratio closer to unity. 

We may assume that if the flexure modes were completely suppressed, the resulting 
frequency spectrum could be approximated by drawing nearly horizontal curves through 
the inflection points of the curves of Fig. 2. This approximation to the frequency ratio 
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Fic. 2. Frequency ratio as a function of length-to-thickness ratio for different degrees of contouring of 
AT-cut crystal. 


w/@ (sometimes called the “frequency constant’’) is often made, in order to show the 
gross behavior of the thickness shear frequency as a function of //h. A set of such curves 
for the five values of r is shown on Fig. 3. This shows much more clearly how the com- 
puted effect of increased contouring is to bring the frequency ratio close to unity. 

The effect of contouring on the first even overtone mode is shown in Fig. 4. This is 
the first overtone mode for which u is an even function of x. Referring to Fig. 1, the 
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Fic. 3. Approximate frequency ratio for different degrees of contouring of AT-cut crystal. 


motion of the surfaces is still parallel to the xz-plane, and the displacement w is still 
proportional to y. However, in this mode the center of the top surface moves in opposition 
to its edges; the displacement wu is a function like y cos (2rz/l). In Fig. 4 the flattened 
portions of the curve with a frequency ratio near 1.018 correspond to this mode. It is 
seen that for this overtone mode the computed effect of contouring is less than for the 
fundamental. Thus the results show that the separation in frequency between the 
fundamental and the first even overtone mode is slightly larger for a contoured quartz 


plate. 
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Fig. 4. Frequency ratio for two degrees of contouring showing the first even overtone mode. 
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Some general remarks may be made about this solution. First, it is in agreement 
with experiment in that the calculated displacement u on the surface does decrease as 
one proceeds from the center towards the edges. Though at any constant height y the 
displacement is greater near the edge, the decreasing height of the surface as one moves 
away from the center is more than enough to offset this. Second, the approximation 
(7) makes possible the simple solution (10) of the differential equations. The relation 
(10) essentially says that the fundamental thickness shear mode, for which u is an even 
function of z, couples only with alternate flexural modes for which v is an odd function 
of x. This simple relationship does not hold when the approximation (7) is not made. 
Nonetheless, experimental evidence seems to show only this coupling with alternate 
flexural modes. 
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BOOK REVIEWS 


(Continued from p. 154) 


The first chapter discusses the “general decision problem”’ a viewpoint that appears to have become 
canonical in most statistical discussions today. Various notions of optimum procedures, Bayes and 
minimax procedures as well as the concepts of invariance and unbiasedness are introduced. Chapter 2 
sketches out the appropriate background in probability theory required for the book. The notion of 
sufficiency is discussed in a measure theoretic background. From this point on the following chapters 
discuss and elaborate in some detail the approaches that have been indicated in the first chapter. Uni- 
formly most powerful tests are discussed in chapter 3. The Neyman-Pearson lemma and extensions 
thereof are introduced. A variety of contexts in which uniformly most powerful tests exist are noted. 
Unbiased tests are examined in the next two chapters. Invariance and a discussion of linear hypotheses 
follow in chapters 6 and 7 respectively. The minimax principle is examined in chapter 8 and the book 
concludes with an appendix. 

The book should be of considerable interest to those interested in a formal and detailed exposition 
of theory of testing and the criteria that have been introduced in this domain. The author is exceedingly 
well qualified since much of his research has been intensively directed toward the theory of testing. 
The orientation of the book is such that it presents techniques directed towards proving the optimality 
or appropriateness of tests proposed some time ago (perhaps partly on intuitive grounds) rather than 


as tools for work in new areas. 
M. RosENBLATTI 


Statistical independence in probability, analysis and number theory. By Mark Ixac. Pub- 
lished by The Mathematical Association of America, distributed by John Wiley & 
Sons, Inc., New York, 1959. xiv + 93 pp. $3.00. 


This little volume is based on the Hedrick lectures given by the author during the summer of 1955. 
The reader is assumed to have some familiarity with the Lebesgue integral and elements of Fourier 
analysis and number theory. The object of the author is to show how probabilistic ideas, in particular 
that of independence, arise in concrete problems that one would ordinarily think of as outside the range 
of probability theory. The first chapter introduces the Rademacher functions as an example of inde- 
pendent functions as well as the usual statistical model of a fair coin-tossing game. In chapter two, 
Borel’s law on the normality of real numbers and convergence (or divergence) of series with random 
signs are considered. The central limit theorem is briefly discussed in chapter 3. Probabilistic ideas 
are applied to some number theoretic problems in chapter 4, which is called “Primes play a game of 
change.’’ The concluding chapter discusses the ergodic theorem as suggested by kinetic theory and 
then applies the ergodic theorem to a problem on continued fractions. The author succeeds admirably 
since the book gives the reader a stimulating contact with the ideas of probability theory and their 
power in investigating problems in other fields. There is a strong emphasis on the analytic aspect of 


probability theory rather than measure theoretics. 
M. RosENBLATT 


Confluent hypergeometric functions. By L. J. Slater. Cambridge University Press, New 

York, 1960. ix + 247 pp. $12.50. 

This book contains a detailed and careful discussion of all properties of the confluent hypergeo- 
metric functions, including Kummer and Whittaker functions, which a mathematical physicist is 
likely to require. The author describes the basic differential equations satisfied by the functions, derives 
the latter’s differential and integral properties and outlines the underlying functional analysis necessary 
in the development of the asymptotic expansions. She presents related differential equations and special 
eases of the functions, such as the Coulomb and Schrédinger wave equations, and concludes with a 
table of 120 pages, calculated on the EDSAC I, of the function ,F[a; b; z] over those ranges most likely 


to be useful to the compute! There is an extensive bibliography. 
WALTER F REIBERGER 
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THE KRON METHOD OF TEARING AND THE DUAL 
METHOD OF IDENTIFICATION* 


BY 
A. I. WEINZWEIG 
University of California, Berkeley 


1. Introduction. Over the past several years, Gabriel Kron has published a series 
of papers expounding his method of solving network problems by tearing the network 
into smaller components, solving the problem on each component, then interconnecting 
the solutions to obtain a solution to the original problem [3, 4, 5, 6, 7]. We wish to present 
a precise mathematical formulation of this procedure. This not only establishes the 
validity of the method but simplifies and extends it, and moreover, leads to a dual 
method we call the method of identification. 

We first formulate a general network problem and establish a necessary and sufficient 
condition for the existence of a unique solution. This has independent interest for it 
simplifies and extends the Kron-LeCorbeiller mixed method of solving network problems 
(2, 8]. Following Weyl and Eckmann [1] an electrical network is considered as a 1- 
dimensional cell complex and the problem formulated in terms of the chains and cochains 
of this complex. The solution of the problem is essentially effected by inverting a certain 
matrix, the matrix of the solution. The method of tearing (identification) transfers the 
problem to a second network obtained by tearing (making identifications in) the original 
network. There the solution matrix is inverted by inverting two matrices, the com- 
ponent matrix and the connection matrix (“interconnecting the solutions”). Although 
the rank of the component matrix is greater than that of the solution matrix, it is strongly 
diagonal and can be inverted by inverting each of the diagonal submatrices (“solving 
the problem on each component”). This is actually a special application of a more 
general procedure developed in Sec. 8 whereby the solution matrix is inverted by invert- 
ing two other matrices, the first of rank greater and the second of rank less than the 
solution matrix. If the inverse of the first is known or for some reason more easily com- 
putable (as in the case of tearing and identification) then this leads to a simpler solution. 
This also furthers Kron’s goal of “storing solutions’. 

For a detailed discussion of the history of the network problem the reader is referred 
to Roth [10, 12]. For an evaluation of Kron’s method of tearing we again refer to Roth 
[12] in addition to the papers of Kron. 

2. The network equations. An electrical network K can be considered as a 1-di- 
mensional cell complex. K is assumed to have lumped design constants with no im- 
pedanceless or admittanceless branches. The k-dimensional chains of K with coefficients 
in the field of complex numbers C,(K) is then a vector space and the k-cochains C*(K) 
with the same coefficients is the dual space. The boundary operator @ is a linear trans- 
formation with dual the coboundary operator 6. Orienting A, the positively oriented 
k-cells can be regarded both as elements of C,(K) and C*(K) and as such define dual 
primitive bases. 


*Received September 29, 1958; revised manuscript received October 9, 1959. 
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The current flowing in each branch in the direction of orientation defines a 1-chain 
2 of K. In the same way, the emf in each branch, the voltage drop across the passive 
coil in each branch and the potential difference across each branch define 1-cochains e 
V and E respectively. To complete the picture, the current flowing out of each node 
or vertex and the potential of each vertex define respectively a 0-chain J and a 0-cochain 
P. These quantities are all illustrated on a representative branch ab oriented from a 
to b (Fig. 1). Clearly V = E+. 
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The impedance matrix [8, p. 2] Z = [Z;,;] and the admittance matrix Y = [}"'’’] represent, 
in terms of the primitive bases, linear transformations ¢:C,(K) — C’(K) and 7: C'(K) > 


C,(K) respectively, where under the usual conditions ¢~ = 7. Kirkhoff’s current and 
voltage laws are expressed by di = J, 6P = E and Ohm’s law by ¢() = V or n(V) = 7. 
These together with V = E + e constitute the network equations. 

3. Canonical decompositions. Let A, be any subspace of the vector space A. There 
is a (not unique) complement A, of A; in A such that A = A, + A, is a direct sum 
decomposition. This decomposition induces the dual decomposition of the dual space 
A* = A* + A*, where A‘ (A‘*), the annihilator in A* of A,(A,) may be regarded as 


the dual space of A,(A,). We denote by o(A,): A, — A the inclusion of A, in A and 
4(A,): A — A, the projection of A onto A, . Then o(A,)*: A* — A* is the projection of 
A* onto A* , that is o(A,)* = 2(A*). 

The linear transformation §: A — A* 
a(A*)to(A,) is non-singular for every subspace A, of A. This implies that ¢é is itself 
non-singular and that £' is also inherently non-singular where we identify A and A** 
under the usual canonical isomorphism. 

Let Z,(K) and B,(K) be the null space and image space of 0. Any decompositions 
C,(K) = Z,(K) + R,(K) or C.(K) = D,(K) + B,(K) are said to be canonical. In this 
case, the restriction of 0 to R,(K) defines an isomorphism of 2,(K) onto B,(A). Observe 
that Z,(K) and B,(K) are canonical. If C'(K) = D'(K) + B'(K) and C°(K) = Z°(K) + 
R°(K) are the corresponding dual decompositions then Z°(K) and B’(K) are precisely 
the null space and image space of 6. 

4. The network problem. The data for a general network problem consists of 
(a) a network K; (b) decompositions C;(K) = C; + C., Co(K) = D, + D, consistent 
with 0, that is 9C; C D; ,i = 1, 2; (c) ip eC. , e, e C¥ , J, e D, , P2 e D§ where C'(K) = 
c* + C* , C°(K) = D* + D* are dual to those given above; (d) a linear transformation 


is said to be inherently non-singular if 
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¢: C,\(K) — C'(K) (the impedance form) or : C’'(K) — C,(K) (the admittance form). 
A solution to the network problem consists of 7, J, EF, e, P satisfying the network equa- 
tions di = 7, 6P = E, n(E + e) = itorg(t) = E +e, such that r(C.)i = 2, , r(C¥) e = 
e,, a(D,)I= I, , «(D%)P = P,. 

Observe that £ or » are dependent only on the design constants of the network and 
not on the network itself. Thus, ¢ or 7 will be the same for any network made up of the 
same coils. 


Theorem. A necessary and sufficient condition for the existence of a unique solution 
to any network problem is that ¢ and hence 7 = ¢“" be inherently non-singular. 


Proof. It follows from the consistency of the given decompositions with 0 that we 
can find decompositions C,(K) = Zy, + Ru + Zi2 + Riz, Co(K) = Dan + Bu + 
Doz + Boo such that C, = Z,, + Ri; , D; = Do, + Boy, it = 1,2 and C,(K) = (Zi: + 
Zi.) + (Ri + Ruy), Co(K) = (Dar + Doz) + (Bo: + Boz) are canonical decompositions. 
Let C'(K) = D" + B" + D”® + B” and C°(K) = Z" +R" + Z” + R” be the corre- 
sponding dual decompositions. In terms of these decompositions the pertinent chains 
and cochains have representations asi = j, + J, +7, + J2,7 =1,+12,P = Pit 
P,, E = E, + E, and e = e, + @i2 + @2 + Coo , Where t2 = jo + do , 
I,e B,;, E; e BY’, P; 2 D;,7 = 1, 2, e = e1 + e2. and E, = 6P,. Then J, = o&, and 
since the restriction of 8 to R,, + Riz = R,(K) is an isomorphism of R,(K) onto B,(K) = 
B. + Ba, J, = OI — dy. 

Setting A, = Z,,, As = Ru + Zi. + Riz yields a decomposition C,(K) = A, + Az 
wherein the component of 7 in A, , namely 7, + J, and the component e,, of V in A¥ = 
D" are known. This is called the solution decomposition. The network problem is solved 


once either 7 or V is known for then ¢(4) = V or V = n(¢) gives the other 
and e = 7(C#)V — E, +e, E = E, + x(Ct)V —e,, P = &'(E) and! = dtisa 
solution. Hence the problem reduces to determining j, or V, = V — e,, . In the ad- 
mittance form V, = [r(A,)no(A*)]"' (J; + t2 — w(Az)n(e:;)) and in the admittance 
form j, = [r(A*)&o(A,)]"" (e:: — 2(Ai)&(é2 + J,)). The transformations (A,) = 
w(A,)no(A*), €(A,) = 2(A*)fo(A,), the solution transformations are non-singular since 


¢ and hence 7 are inherently non-singular. 

The uniqueness follows readily, for any other solution 2’, E’, I’, P’, e’ can be repre- 
sented in terms of the decompositions above as 7’ = jf + Jf + jo + J2,¢€ =e, 4+ 
ero te ten, E=E+E,P =Pi+P,,1 = 1, + If whence If = dJ, = 1,, 
Ji = a ') — J, = J, , ete. We remark that P is uniquely determined only to within 
a constant value on each connected component of K. Knowing the potential of one 
vertex of each component, say it is grounded, P can be uniquely determined. 

To establish the necessity of the condition, let C,; be any subspace of C,(A) and 
i, e C, such that r(C*)fo(C,)z, = 0. Passing to a new network K’ if necessary, C, can 
be taken as a subspace of Z,(K). Such a network can always be obtained b:  ‘entifying 
vertices in K; (if all the vertices are identified, then Z,(K') = C,(K’)). The> .et R,(K) 
be a complement of Z,(K) in C,(K) and Z,, a complement of C, in Z,(K). Set C, = Z,, + 
R,(K). Then C\(K) = C, + C. , Co(K) = D.(K) + B,(K) are decompositions con- 
sistent with 0 and together with 7, = 0, e, = 0, 7, = 0, P. = O and £, define a network 
problem with solution 7 = 0,e = E = 0, P = 0,J = 0. By the unicity of the solution 
it follows that 7; = 0 so that r(C*)to(C,) is non-singular and ¢ is inherently non-singular. 

The solution outlined in the proof above is actually a generalization and simplifi- 
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cation of the Kron-LeCorbeiller orthogonal or mixed method [2, 8]. Roth [10] con- 
sidered a special case of the network problem, where the solution decomposition is 
‘anonical and established a necessary condition for the existence of a unique solution. 
In a subsequent paper [12] he established a necessary and sufficient condition, namely 
that ¢(z) ¢ = 0 only if z = 0. If ¢ satisfies this condition then ¢ is said to be ohmic. Actually, 
¢ is ohmic if and only if ¢ is inherently non-singular. Let C, be any subspace of( C,K) 
and ze C, . Then o(C,)*fo(C,)¢ = O only if ¢(2)¢ = O for the annihilator of C, can be 
taken as the complement of C* in C’(K). Hence if ¢ is ohmic, 7 = 0 and ¢ is inherently 
non-singular conversely, let ¢ be inherently non-singular and 0 ¥ 7 e C,(K). Then taking 
C,, as the subspace of C, (A) generated by 7, it follows from the fact that o(C,)*fo(C,)i # 0 
that ¢(7)¢ + 0 and ¢ is ohmic. 

5. A more general formulation. The network problem can be formulated more 
generally. Let X and Y be two vector spaces, §: Y — Y a linear transformation. Then 
the roles of C,(K), Co(A) and @ in the network problem as formulated above, are re- 
placed by X, Y and & respectively. Thus the data for a network problem consists of 
spaces X, Y; linear transformations §: XY — Y, \: X — X*; decompositions X = X, + 
X,, Y = Y, + FY, consistent with ¢; and xt e X* , 2,2 X.,y, e Y, , y$ e Y* , where 
X* = X* + X*, Y* = Y* + Y* are the corresponding dual decompositions. A solution 
then consists of x e X, «* « X*, ye Y, y* « Y* such that tx = y, Ax = a* + £*y* and the 
projections of x, 2*, y, y* in X*, X,, Y,, YE are xt, x,y, , y% respectively. Thus we 
have a problem in linear algebra. 

6. Matrix formulation. Choosing bases in (,( 4), C'(K), C,( K), C°( RK), the trans- 
formations 7, {, 0, 6 can then be represented by matrices [7], [¢], [0], [6], the chains 27, 
I and the cochains e, £, P by column matrices [7], [7], [e], [£], [P] respectively, where 
[n]"" = [¢], since ¢ is assumed inherently non-singular. The network equations then 
become the matrix equations [e] + [£] = [¢][z], [7] = [n]({e] + [£)), [e][2] = [7], [é[P] = 
[EZ]. When the bases in C,(K), C'(K) and C,(K), C°(K) are dual then [d]* = [6]. If the 
basis for C,(K) contains a basis for Z,(A), the remaining elements span a complement 
R,(K) of Z,(K) so that C\(K) = Z,(K) + R,(K) isa canonical decomposition. The basis is 
then called a canonical basis for C,(K). Similarly, we define a canonical basis for C,(A). 
The isomorphism @ | R,(K) of R,(K) and B,(K) can be represented by a non-singular 


matrix [d]’. In this case [0] has the form 
_ Fl0) fay’) 


[0] [0] 5 


[d] 


where [0] indicates the appropriate zero matrix. Finally, if the canonical basis for C)(A) 
is such that B,(K) is spanned by the boundaries of those elements of the canonical 
bases for C,(K) which span R,(K), then [d]’ is just the unit matrix. We say the canonical 
bases are compatible. With compatible canonical bases, 


“7 fi) r I), ’ 0 ; ? P 7 
a=)" m=("* w=] Oo} ele), 
[J] LO. [E]._| [P].- 
where [J] = [/], , [E], = [P], . In this way, all the information is transferred directly 
to C,(K) and C'(K) which leads to the solution decomposition C,(K) = A, + A.. 


In practice the data for a network problem is given with respect to some bases. 
The design constants determine both the admittance and impedance matrix which are 
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the matrix representations of 7 and ¢ with respect to dual primitive bases. Let < 6, , +++, 
by > denote the primitive basis in terms of which (and its dual) [n] = Y and let < 67, 

b’ > be the canonical basis which effects the solution decomposition C,(A) = 
A, + A, , that is {b, , --- , b,} spans A, and {b,,, , --- , by} spans A. . The change 
from coordinates with respect to the primitive basis to coordinates with respect to the 
canonical basis is given by a non-singular matrix T = [t,;], 7, 7 = 1, --: , N. Thus if 
[7] and [{2]’ are the matrix representations of 7 in terms of the primitive and canonical 
bases respectively, then [7]/ = 7{z]. Similarly 7*, the transpose of 7' represents the 
change in coordinates in C'(K) from the dual canonical basis to the dual primitive basis. 
Let 7, be the submatrix of 7 defined by [¢;;], 7 = 1, ---,N,j=rt.1-:+,N.T, 
can be regarded as the matrix representation of the projection (A,) in terms of the 


primitive basis for C,(A) and the basis < b/., , --- , bx > for A, . Similarly 7'§ represents 
the inclusion o(A* r(A,)* in terms of the dual bases. The matrix representation 
for the solution transformation n(A,) in terms of the basis < b/,, , +--+ , by > for A» 


and the dual basis for A* is then 7,Y7* . This then is the matrix to be inverted in the 
solution of the network problem. 

The computations can often be simplified by the judicious choice of bases leading 
to simpler matrix representations. This has already been indicated above. In particular, 
the l-simplices of a maximal tree in K generate a complement F,(K) of Z,(AK). Each of 
the remaining simplices gives rise to a simple circuit and the simple circuits so formed 


are a basis for Z,(A). The basis for C,() so constructed is called a simple canonical 
basis. The coordinate transformation 7 from a primitive to a simple canonical basis 
las the form 
o [1] [0] | 
y {Ul OF) 
_ | 
Lv, Tey 
1 


where [1] is the unit matrix and the only non-zero entries of 7, are + 1. Moreover T 
is obtained from 7 by replacing the submatric 7, by — T; . 

lor a more detailed discussion of this the reader is referred to a forthcoming work 
of the author [13]. 

Since the solution of a network problem reduces essentially to inverting the solution 
matrix, we often speak of the inverse of the solution matrix as the solution. 

7. The transformation matrices. For any two networks A, A’ made up of the 
same coils, C;(K) and C,(K’) are canonically isomorphic and may be identified. How- 
ever, to each of the networks K, K’ corresponds a subspace Z,(A), Z,(K’) which is de- 
termined by the topology of the network. 

Consider the set of canonical bases arising from the various networks made up of 
the same NV coils. Observe that a canonical basis for one network need not be a canonical 
basis for a different network. Let b, be any such basis and Y", Z, , 2", e, , 2, the matrix 
representations of », ¢, 7, e, E in terms of 6, and its dual. If 6,, is another canonical basis 
the change from coordinates with respect to b,, to those with respect to b, is given by 
the non-singular transformation matrix C™ . Then, 


“Ce, ann; B.-~ CR, £408, = Cire. 


Moreover C” (C")~' and if bp is another canonical basis then C? = C?C™. The set of 
transformation matrices [7] is a subset of the full linear group of N X N non-singular 
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matrices but they do not, as Kron asserts [7] form a subgroup since the product CC} 
is defined as a transformation matrix only when n = gq. 

If b, and b,, are bases related to the networks K, and K,, , the transformation matrix 
C™ permits the “transference” of the data of a network problem on K, to one on K,, . If 
the problem has there been solved, say 7” and V,, have been determined so that 7” = 
Y"V,,, then” = C™%”", V, = Cz* V,, is such that 


Y°V, = CaY"C.C; V. = CaY"V.. = Cu” = 3% 


is a solution on K, . Thus problems and solutions can be transferred from one network 
to the other. 

8. The use of known solutions. Let C,\(K) = A; + Az be the solution decomposi- 
tion for a given network problem and 7(A,) = 2(A2)no(A¥%), §(A:) = 2(A%) fo(. A,) the 
the admittance and impedance solution transformations respectively. If 7 = ¢ is 
known, the admittance form of the network “egy can be transformed into the imped- 
ance form and can thus be solved by inverting ¢(A,) rather than 7(A,). Whenever the 
rank of the former is less than that of the latter, this will be a simpler problem. In fact, 
we can write 7(A.)~'i, = ¢(i2 — €(A,)~* (w(A*)gi2)). Thus knowing the inverse of 7, 
enabled us to invert n(A,) by inverting instead ¢(A,) 

More generally, let C;(K) = B, + B, with B, D A, and let n(B,) = x(B,) no(B%): 
Bt — B, where C'(K) = B* + B* is the dual decomposition. We can find a complement 
A; of A, in B, so that B, = A; + A, . Then, as above, n(A2)~'(i2) = n(B,)~*(22 
(3(A*)n(B,)~'o(Az))~ 2(A*)n(B,)~'t2) so that, knowing 7(B,)~", n(A,)~' can be com- 
puted by essentially inverting m(A*)(B,)~'o(A;), the connection transformation. (see 


Sec. 9 below). 
1 


There is an analogous impedance formulation. Thus ¢(A,)7'V, = (V,; — (A, 
(x(A*)nV,)) and more generally, let C,(K) = Bi + Bi with BJ D A, , ¢(Bf) = w(Bi*) 
to(B{), where C’(K) = B{* + Bj* the dual decomposition. Then as before, Bf can be, 
represented as BJ = A, + Aj and 

¢(A,)7V, = (BY) (V, — (( A) ¢(B))"0(Af*)) '( A) ¢(B)) VV). 

To obtain a matrix formulation of the above procedure let b, = < bi, --- ,b% > be 
a canonical basis effecting the solution decomposition A, + A, , that is {b} , On 
span A, and {b,*', --- , bX} span A, . In terms of this basis and its dual. 7 is represented 
by the matrix Y”. Partitioning Y”’, 

oe ee cd 

y* = : | 

LY; Yj 
where Yj is an N — r X N — r submatrix of Y" ge ee the solution transfor- 
mation mt 4.). On the other hand, let b,, = < b,, --: by > effect the decomposition 
C,(K) = B, + B, where {b1 , --- , b5} span B, and {b%"', --- , bX} span B, , and let 


Y” represent 7 in terms of the basis b,, and its dual te Y™ be ‘de N-sXWN 

submatrix of Y” which represents 7(B,) in terms of the basis ‘< bi*) , --- , b* > and its 
dual. The basis b, is chosen to effect both decompositions, that is < b%*', --- , bX > 
spans B, and D” the n represents the change in coordinates in B, . In practice, the com- 
putations are carried out in terms of the basis b, while Y% and its inverse are given in 


1 


terms of the basis b,, hence D™’ Y””'D™ gives the matrix representation of »(B, in 
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terms of b, . Finally, if D™ = [d,,;],7,j = s+ 1, ---,N then EZ” = [d,,;],i = s + 1,---, 
r,j = s+ 1, --- , N gives the matrix representation of o(A;) in terms of the bases 
< bf*', --+ , bl > and < b%"", --- , b§ > so that the connection matrix to be inverted 
is E"Y"E*.. 

9. The method of tearing and the method of identification. In both these methods 
the topology of the network is utilized to obtain a decomposition C,(K) = B, + B, with 
n(B,) and hence Y™ above, strongly diagonal. In this way the matrix representation of 
n(Az), Y? is inverted by inverting several matrices of smaller rank. 

Let K, be a network or complex and K, , --- , K, subcomplexes of Ky such that 
every 1-cell of K, lies in exactly one of the subcomplexes K,; , i = 1, --- , k, and no coil 
of K; is inductively coupled to a coil of K; , for i ¥ 7. Then 

C,(Ko) aad C,(K,) + C,\(K2) = Maite - C,(K,). 
Moreover, the admittance transformation 7 respects this decomposition, that is 7, = 
m(C,(K,)) no(C'(K,)) is just the admittance transformation for the network K,; , i = 
1, --- , k. Choosing a (canonical or simple canonical) basis consistent with this de- 
composition, the matrix representation of 7 is strongly diagonal. 

Let K,,, be the complex consisting of the disconnected subcomplexes K, and let 
0,(K,) = Z,(K,) + R\(K,), i = 1, --- , k, be canonical decompositions. Then Z(K,.,) = 
Z,(K,) + --- + Z,(K,) and Z,(K,.,) + R,(K,i.;) is a canonical decomposition for 
(,(K) where R,(K,.,) = R,(K,) + --- + R,(K,). Clearly Z,(K,.;) C Z,(Ko) so that 
R,(Ko) can be chosen such that Z,(K,) + R,(Ko) is a canonical decomposition for 
C,(Ko) = C,(Ky.,) and R,(Ko) C R,(K,.,). Setting »(K;) = 2(R,(K;) no (B'(K,)), 


n(K,.,) agrees with n(K;) on B'(K;) for i = 1, --- , k. Hence for the proper choice of 
basis, the matrix representation of n(K,.,) is strongly diagonal. Thus, »(K,.,) can be 
inverted by inverting each of n(K;), i = 1, --- , k. If the subcomplexes K; ,i = 1, --- , & 


are chosen to consist of as few distinct complexes as possible, the matrix representations 
of »(K,) for the “same’’ complexes, can be made identical thus simplifying the task of 
inverting n(K;,.,;). 

If the solution decomposition C,(Ky) = A, + A; for a general network problem on 
K, is such that A, C R(K,.,) then the procedure of Sec. 8 can be applied. This is the 
method of tearing. If dim A, = N — k and dim R,(K,.,) = N — s, then the inversion 
of the N — r X N — r solution matrix is effected by inverting each of the diagonal 
submatrices of the strongly diagonal N — s X N — s matrix and ther — s Xr—s 
connection matrix. 

The problems considered by Kron illustrating his method are all in admittance form 
with canonical solution decompositions (and hence of the special form considered by 
Roth [10]). In the method of tearing, the problem is transferred from Ky, to A,,, where it 
is now of the more general form treated above, for the canonical decomposition for Ko 
is not a canonical decomposition for K,., . In this case, however, A, = R,(K.) C Ri(Ay.1) 
and the above considerations apply. In general, however, A, D R,(Ko) so that the 
condition that R,(K,.,;) D A2 may not be met. The problem can then be cast into the 
impedance form by inverting the strongly diagonal matrix representation of 7. Then 
A, C Z,(K,.,) and we can apply the admittance form of the technique of Sec. 8 wherein 
we invert the solution matrix by inverting the strongly diagonal matrix representation 
of €(Kys;) = 2(D'(Ky.,)) to(Z,(Ky.,)) and an impedance connection matrix. 

The greater the number of pieces into which K, is torn to form A,., , the smaller 
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dim Z,(K,.,) and the greater dim R,(A,.,), and hence the rank of the connection 
matrix. The smaller the individual pieces the smaller the rank of the transformations 


nik 
the technique of Sec. 8, R,(A,.,) must contain A, . 


T 


In practice, these two things have to be balanced. However, in order to apply 


he method of tearing is not in general suited to problems in the impedance form 


for in this case, to apply the technique of Sec. 8 we want to shift the problem from 


nt 
can be chosen strongly diagonal. Since Z,(K») D A, , it is sufficient to have Z,(K,..) D 


ZAK 
of the subcomplexes K; , 7, 7 = 1, 


two 


modified b: 


obtained 


Ne 


Ba 
1 
9 ( 
S, 
4.G 
5 
6. ¢ 
7. 
8. | 
9. J 
10. J 
ye 
12. J 
13 


Vi 


»f 


a new network A,.. where Z,(K,..) D A, and the matrix representation of ¢(K,.>;) 
K,.2 is obtained from K, by identifying in A; , all those vertices common to 

= , k. This process of identification can be 
y selective identification provided only that if K‘ is the subcomplex of K;,,. 
by identifications in K, then K,.. can be torn into the disjoint subcomplexes 


thout opening any cycles of A, 


i 
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A CONTINUUM MODEL FOR TWO-DIRECTIONAL TRAFFIC FLOW* 


BY 
J. H. BICK** ano G. F. NEWELLt 


Brown University 


Abstract. The flow of traffic in two directions of an undivided highway is investi- 
gated using the equations of continuity and assumed empirical relations between the 
average velocities and densities in both lanes. These lead to a pair of quasi-linear partial 
differential equations. Even if the velocity in one lane depends only very weakly on the 
density in the other lane, it is found that for a certain small range of densities the equa- 
tions are of elliptic rather than the expected hyperbolic type. For densities outside this 
range, solutions of the equations can be found for various special types of initial con- 
ditions. 

1. Introduction. A number of papers have now been published treating the theory 
of traffic flow from the standpoint of the dynamics of individual cars, statistics and 
continuum mechanics. Lighthill and Whitham [1], Richards [2], and De [3] have studied 
unidirectional flow of traffic by assuming that the empirical relation between the average 
velocity of cars and the density of cars, determined under steady state conditions, 
applies also to non-steady flows. They then use the one dimensional equation of con- 
tinuity to describe how the densities vary with time for various types of disturbances. 

Here we shall attempt to extend their work to include two lanes of traffic flowing 
in opposite directions. We use two equations of continuity and assume that the average 
velocities are functions of the densities in both lanes. If we let p and g denote the densities 
of cars and U and V the average velocities in the right-moving lane (p-lane) and the 
left-moving lane (q-lane), respectively, the equations of continuity for the two lanes are 


p. + (pU), = 0, q. + (qV), = 0, (1) 


in which subscripts denote partial differentiation with respect to the subscript variable. 

We will assume that p, g, U and V satisfy the following conditions: (A) U and V are 
functions of p and q only; i.e., the highway is uniform in both space and time; (B) p and 
q are never negative; (C) U is non-negative and V non-positive; (D) U and V have 
continuous first derivatives with respect to p and q; (E) U is a monotone decreasing 
and V a monotone increasing function of both p and q; (IF) for each value of q, there is a 
value of p at which U = 0 and for each p there is a value of g at which V = 0; (G) 
U(p,q) = —V(q, p);i.e., both lanes are physically identical. Property (F) in conjunction 
with (A) and (B) implies that the allowable domain of solutions is finite and bounded 
by the four curves p = 0, q¢ = 0, U = Oand V = O in the (p, q) plane. 

If shocks develop, they must satisfy the integrated equations of continuity 


A[U(A, B) — W) = p[UQ,9 — W), BIV(A,B) -—W)=aIV~.9—-W)], (2) 


*Received Sept. 15, 1959. 

**Now with Atomics International, Canoga Park, Calif. The work of J. H. B. was supported by 
the International Business Machine Corporation of New York City. 

tAlfred P. Sloan Research Fellow. 








192 J. H. BICK AND G. F. NEWELL [Vol. XVIII, No.2 


where W is the shock velocity, A and B represent the densities in the p- and q-lanes 
respectively on the left side of the shock and p, q the corresponding densities on the right 
side of the shock. It follows from this, that if there is a shock in one lane, there must 
also be a contiguous shock in the other. If, for example, A were equal to p (no shock 
in the p-lane), Eq. (2) gives U(A, B) = U(A, q). By property (£), this implies that 
B = q, so there can be no shock in the qg-lane either. 

Several methods exist for solving quasi-linear differential equations such as Eq. (1). 
First, there is the method of characteristics in the physical or (x, #) plane. The character- 
istics of a system of equations are those curves in the (x, ¢) plane along which disturbances 
will propagate and a set of two first order partial differential equations is called hyperbolic 
if there are two real, one-parameter families of such curves. If, for hyperbolic equations, 
initial conditions are specified along a non-characteristic curve, the solution, can be 
obtained by numerical integration along the characteristics [4, 5]. 

A second method of solution results from exchanging the roles of dependent and 
independent variables. The technique is well known in the theory of compressible fluids 
and is referred to as the hodograph transformation [6]. The equations of continuity in 
the (p, g) or hodograph plane are, 


z, — (U + pU,)t, + pU,t, = 0, (3) 
x, — (V + qV2)t + qV ot, = 0. , 


Since this system is linear, the characteristics, unlike those in the physical plane, are 
independent of initial conditions. The initial conditions, given by specifying p and q 
as functions of x at time zero, describe a parametric representation of a curve, p = 
p(x, 0), g = q(x, 0), in the hodograph plane along which x and ¢ (¢ = 0) are specified. 
Such data, if it does not lie along a hodograph characteristic can be uniquely extended 
to define a solution z = x(q, p), t = t(p, q). 

Although one would expect the equations of continuity to be of hyperbolic type, the 
only type which gives a unique solution for initial data p(x, 0) and q(x, 0), Eq. (1) is 
not always hyperbolic. The physical plane characteristics are defined by the velocities 


Qw =2 = = (pU), + (qV), £{[(pU), — (qV).]) + 4qpU,V,}'”. (4) 
For every fixed value of g, pU is zero when p is zero and also when p has its maximum 
value, i.e., when U is zero. It is positive between these points, and therefore, by virtue 
of property (D), it has a maximum at which (pU), is zero. Similarly, gV has at least 
one negative minimum in the allowable range of q for every fixed value of p. Conse- 
quently, there exist values of p and q for which [(pU), — (qV),]’ vanishes. By properties 
(B), (D), and (E), the term 4qpU,V, is negative for all g and p greater than zero. For 
a certain range of p and q, the characteristics are therefore complex, and the equations 
of continuity are of elliptic type. Outside this region the equations are hyperbolic. 

As yet we have found no satisfactory explanation of why the equations should be 
elliptic nor any satisfactory suggestion as to what one should do about it. Although 
one could imagine hypothetical situations that would give rise to densities in the region 
of ellipticity, it is not obvious whether or not such situations could be realized in practice. 
It may, however, also be true that one is not justified in assuming that the relations 
between the velocities and densities, determined under stationary conditions, are also 
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valid under non-stationary conditions, in which case much of the theory described here 
as well as that in reference [1, 2, 3] may collapse. One might, in fact, seriously question 
any theory in which the future behavior of a system is determined with no reference to 
the distribution of the velocities of individual cars. On the other hand, the ellipticity 




















0 


Fig. 1. The bounded domain of solutions in the hodograph plane and the region of ellipticity obtained 
from Eq. (5a) are shown for 8 = 14. The two families of characteristics are distinguished by the labels, + 
or —, and the arrows indicate the direction of increasing wave velocity. 


arises when the densities are such that, in the absence of a coupling, the wave velocities 
in the two lanes are nearly equal. Under such conditions, one could reasonably expect 
any coupling between the lanes to produce some unusual effects. 

2. A special case. It has been found experimentally for unidirectional traffic that 
the average car velocity is approximated reasonably well by the linear relation 
U = U, — ap. To illustrate some of the consequences of the present model, we consider 
the velocity-density relations 


U = Vo — ap — 84; V = —U, + ag + Bp, (5) 


where a and @ are constants with a > 8 > 0. These satisfy all the properties listed in 
the previous section and reduce to Richards velocity functions if 8 = 0, p = Oorg = 0. 
Furthermore they have the form of the first three terms of a Taylor’s series expansion, 
so if we know the velocities for some values of p and g, we could use Eq. (5) for small 
variations in p and q about these values. 

The velocity relations, the equations of continuity, and the shock equations, may 


be non-dimensionalized by the substitutions U’ = U/U, , V’ = V/U,, p’ = ap/Uo, 
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q = aq/l,, B’ = B/a, t = Ut, and W’ = W/U, . Dropping the primes for the sake 
of simplicity in notation, we obtain from Eq. (5) 

U=1-p- £q, V = -1+q+ Bp, (5a) 
whereas Eqs. (1-4) remain unchanged. 

The region of the hodograph plane in which the equations of continuity are of elliptic 
type is the interior of an ellipse (see Fig. 1) obtained by setting the radical of Eq. (4) 
equal to zero. This ellipse is tangent to the two axes of the hodograph plane and to the 
line p + gq = 1. It has a minor axis of width (2)'*6/(2 + 28). 

The characteristics in the hodograph plane are determined by the differential equation 

Bp (dq) — [2 — (2+ B)(p + q)] dp dq + Bq (dp) = 0 (6) 
the solution of which can be found by means of the Legendre transform [7]. In parametric 


form, the solution is 


~~VWX)_ yy 7a ; 
p=xX AX Y(X), a. (7) 
with 
om i hay Of dy 
V(x) =f —_ “ __ 4G}, 7% 
l+~, 1 xA+D=™ as 
e = B/(2 + 28), (7b) 


and C is a constant of integration. For any X, Eq. (7a) determines Y and Eq. (7) then 
gives both p and q as a function of the parameter X. In addition Eq. (6) has the singular 


solution 

ergs (8) 
which forms an envelope of the general solutions. Some characteristics for 8 = 1/3 
are shown in Fig. 1. A more detailed description of the hodograph characteristics par- 
ticularly for 8 < 1 is given in [8]. 

The shock relations, Eq. (2), consist of two equations in the five variables A, B, p, ¢ 
and W. If we eliminate W from these equations, we obtain a quadratic equation in p, 
q, A, and B which for any fixed values of A and B can easily be solved for g as a function 
of p or vice versa. The resulting curve in the (p, qg) plane is known as the shock polar 
corresponding to the given values of A and B. If (A, B) lies in the hyperbolic region of 
the hodograph plane, the shock polar will pass through the point (A, B) twice and will 
be tangent each time to one or the other of the two hodograph characteristics through 
the point (A, B). Some typical shock polars are shown in Fig. 2. 

Although it is always possible, for some shock velocity W, to find a shock satisfying 
Eq. (2) between any denisities (A, B) and densities (p, g) lying on the appropriate shock 
polar, the resulting shock may not be stable (even for 6 = 0). If at zero time, we specify 
continuous strong discontinuities that lie on a shock polar and we wish to calculate the 
densities at time Af later, Eq. (2) gives us two equations in the five unknowns A, B, 
p,q and W at the time At. Additional relations are given by Eq. (1) applied along the 
physical plane characteristics provided these characteristics emanate from the line 
t = 0 along which the initial data are specified. It follows that if a strong discontinuity 








1960] TWO-DIRECTIONAL TRAFFIC FLOW 195 








a (AB) “te 
| 








. (A,B) 1? 
(a) (b) 


























(c) (d) 
Fic. 2. The shock polars of points, (A, B), above and below the ellipse are drawn for 8 = 0 and < 1 
in (a) and (b). The dashed lines are the shock polars for 8 = 0, and the solid lines, for 8 «< 1. Shock 
polars for a point, (A, B) located just above and just below the line p + g = 1 are shown for 6 = 4 


in (ec) and (d). The stable segments of the shock polars are shown by the heavy lines. 


in initial data is to persist, three of the four characteristics issuing from the line ¢ = 0 
on either side of the shock path must intersect this shock path in the (x, é) plane. If less 
than three characteristics intersect the shock path, we should look for wave type solutions 
of Eq. (1). If four intersect the shock path, the problem is overspecified. Physical plane 
characteristics originating from the left side of the shock will meet the shock path if 
their velocities are greater than the shock velocity W, whereas those on the right side of 
the shock path will meet it if their velocities are less than the shock velocity. It is possible 
to distinguish two types of stable shocks accordingly as two of the three characteristics 
that must cross a stable shock originate from the left or from the right of a shock. We 
refer to these as left or right shocks, respectively. If a shock is to be a stable left shock, 


the inequalities 

w.(A,B) > w(A,B) > W > w(p,q, wilp,g>W (9) 
must be satisfied and to be a stable right shock, the inequalities 

wi(A,B) > W > w.(p,q) > wp, gq), w.(A,B)< W (10) 


must be satisfied where w, and w_ represent the positive and negative roots of Eq. (4). 

Since the stability criterion adopted here is based upon the properties of the physical 
plane characteristics, we can say nothing about the stability of shocks when either of 
the points (A, B) or (p, q) lies in the elliptic region of the hodograph plane. 
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3. Discontinuous initial data. In certain cases, as for example when a stoplight 
changes from red to green, we may expect the initial conditions to be discontinuous. 
We shall consider here solutions of the equations of continuity for initial data of the 
type p(z, 0) = p: , g(x, 0) = gq, for x < O and p(z, 0) = p, , g(x, 0) = q, forz > 0 
where p; , ?, , 9: , and g, are constants. Such initial data map into the hodograph plane 
as just two points, (p; , g:) and (p, , q,). 

The equations of continuity and these initial conditions are invariant under the 
continuous group transformation x’ = yz, t/ = yt, with y any arbitrary constant and so 
we may expect to find solutions in which the dependent variables, p and gq, are functions 
of z/t only [9]. In this case, an increase in time corresponds only to a rescaling of the 
x axis. Thus, if a solution exists, its form will be the same for any positive time. Sub- 
stitution of y = x/t into Eq. (1) gives, 


dp | (11) 


ee a 
G=-g~-—- Se; G+e~— 39 BP) ay a 


dy dy 
Since Eq. (11) is homogeneous and linear in the derivatives, either dp/dy = dq/dy = 0, 
in which case the solutions, p and q, are constants or the determinant of coefficients must 


vanish i.e., 
y = (1 — 8/2)(q-—p) + {11 —(1 + 8/2)(p + OP — Bpgq}'”. (12) 


In the latter case, substitution of Eq. (12) into Eq. (11) gives the differential equation 
for the hodograph plane characteristics, Eq. (6), which has been solved in Sec. 2. Thus, 
for each hodograph characteristic there exists a special solution p(y), g(y) obtained 
by assigning to each point (p, g) along this characteristic the value of y given by Eq. 
(12). Whether one chooses the + or — sign in Eq. (12) depends upon which of the two 
families of characteristics the one in question belongs. The characteristics in Fig. 1 are 
labeled + or — accordingly as one should use the + or — sign in Eq. (12) and we shall 
refer to the corresponding solutions of Eq. (11) as plus or minus waves, respectively. 

The complete solution of Eq. (11) is obtained now by considering all possible ways 
in which one can decompose the y-axis into intervals on each of which p(y) and q(y) 
are either constants or form plus or minus waves and at the boundaries of which p(y) 
and q(y) are either continuous or have stable shocks. For any shock or wave, y is finite, 
so there must exist a finite y to the left of which p and q are constants, p = p, and g = q:, 
and another y to the right of which p = p, and gq = q,. 

The values of y given by Eq. (12) are the velocities of the physical plane character- 
istics corresponding to given values of p and q, as one can see by comparing Eq. (12) 
with Eq. (4) and (5a). The conditions for stability equations (9) and (10), therefore, 
take the special form in which w., (p, q) is replaced by y = y, (p, g) evaluated from Eq. 
(12) with the positive root, etc., and W is the value of y for the shock. 

There are only a few combinations of waves and/or shocks that can be used to con- 
struct a permissible solution without violating the shock stability conditions or the 
requirement that p and q be single valued. Suppose, for example, that as y increases, a plus 
wave were to terminate at some point where the densities are (A, B) and the velocity y is 
y, (A, B). One could attach to the right of the plus wave, y > y, (A, B), a region of 
constant densities but according to Eqs. (9) and (10), a shock from the densities (A, B) 
to any other densities (p, g) would have a velocity W which is less than y, (A, B), i.e., 
it would occur to the left of y, (A, B). A minus wave starting at a point where the 
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densities are (A, B) would also have a velocity y_ (A, B) < y, (A, B), to the left, of 
y. (A, B), unless per chance (A, B) lies on the ellipse where y_ = y, . If (A, B) does lie 
on the ellipse, there is the possibility that a plus wave can be followed by a minus wave, 
provided that the plus wave characteristic through (A, B) approaches (A, B) with y 
increasing and the minus wave leaves (A, B) with y increasing. Figure 1 shows the 
direction of increasing y for the various characteristics and one can see that the only 
characteristic having this property is the singular characteristic, Eq. (8), and the only 
point (A, B) where this can happen is the point (1/2, 1/2). We conclude from this that 
the only thing that can appear to the right of a plus wave is a region of constant densities 
unless the plus wave belongs to the singular characteristic and terminates with densities 
(1/2, 1/2). 

If one has a stable right shock from any density (A, B) to densities (p, g), it can be 
followed by nothing except a region of constant density because, by Eqs. (9) and (10) 
again, the shock velocity is larger than either y, (p, q) or y- (p, q). Similarly one can 
show that from a left shock, one may go to a region of constant densities which in turn 
may be followed by either a plus wave or a right shock but not another minus wave 
or left shock. 

It follows from this that the only combinations of waves and shocks that can form 
permissible solutions are (1) a single wave or shock going from (p; , g:) to (p, , g,) with 
regions of constant density on either side, (2) a minus wave of left shock going from 
(p: , 9x) to some densities (p, , g;) followed by a region of constant densities (p,; , g;) and 
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Fia. 3. Points in numbered regions, shown for 8 = 0.1, can be reached by various allowed solution 
sequences starting from some point 0 below the ellipse. The dotted curves represent the maps of sample 
solutions. 
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then either a plus wave or a right shock to (p, , g,) (3) a minus wave or left shock going 
from (p; , g,) to a point on the singular characteristic followed by a region of constant 
density, then the plus wave to minus wave sequence described above, another region 
of constant densities, and then possibly also a right shock or plus wave to (p. , q,). 

Figure 3 shows the various domains in the hodograph plane that can be reached by 
these allowable combinations of waves and/or shocks from a point (p, , q,) designated 
in Fig. 3 as the point 0 located below the ellipse. The line 0A is the image of a minus 
wave, a hodograph characteristic through 0 in the direction of increasing y—, and the 
line OC represents a plus wave, a hodograph characteristic through O in the direction 
of increasing y, . The lines OD and EK depict those segments of the shock polar which 
are the loci of stable left shocks from the point 0, and the lines 0/ and L.A/ represent 
stable right shocks from 0. Points on these lines can be reached by a single wave or shock 
from 0. 

Any point in region | of Fig. 3 can be reached from 0 by a combination of a minus 
wave along 0A followed by a plus wave into region 1. The boundary AB of this region 
is the plus characteristic through A. A sample path is sketched in Fig. 3, and the physical 


in Fig. 


plane solution corresponding to it is shown in Fig. 4(a). All points in regions 2: 
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Fic. 4. The two solutions whose hodograph maps are shown in Fig. 3 are represented by plott ind q 
as functions of y = x/t. Figure 4 (a) shows a minus wave followed by a plus wave and Fig. 4 i left 
shock followed by a plus wave. 

3 can be reached by a combination of a left shock along the line OD followed by a plus 


wave into region 2a. A sample path is shown in Fig. 3, and the corresponding physical 
plane solution in Fig. 4(b). All points in the region 2b can be reached with a left shock 
from 0 to the line EK followed by a plus wave into region 2b and all points in region 3 
can be reached with a minus wave followed by a right shock. 

Points in region 4a can be reached with a left shock along the line OD followed by 
a right shock, and points in 4b by a left shock to the line EK followed by a right shock. 
The line WK is the locus of the end points of the stable shock polars for right shocks 
from OD. As one approaches any point on the boundary KF between these regions, one 
finds that the velocities of the two shocks used to reach this point become equal and the 
two shocks become a single shock. The line KF is, therefore, part of the shock polar 
determined by 0 but it is an unstable part because four physical plane characteristics 
cross the shock path rather than the usual three. 

If we set initial data with (p, , q:) at the point 0 as in Fig. 3, we may set (p, , q 
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anywhere in regions 1 to 4 and still obtain a stable solution. For any given (p, , q:), 
however, stable solutions are not assured for all values of (p, , g,). Regions 1 to 4 do 
not cover the elliptic region or even all of the hyperbolic regions of the hodograph plane. 

When there is no coupling between lanes (8 = 0) the densities are independent of each 
other, and every point in the hodograph plane may be reached from every other point 
by the solutions of two Richards-type problems. That there exist regions in the hodo- 
graph plane that cannot be reached from an arbitrary point (p, , g:) by any allowable 
combination of two shocks and or waves for 6 > 0 is due to the fact that these particular 
Waves cannot be extended through the ellipse. 

The diagrams analogous to Fig. 3 for various other positions of the point 0, such 
as above the ellipse or between the ellipse and the singular characteristic, will have 
quite different forms. Whereas there was no opportunity to use the singular character- 
istic in any solution sequence of Fig. 3, it does appear for some other positions of the 
point 0. That one can not reach every point in the hodograph plane with an allowable 
sequence of waves and shocks from 0 is, however, true for all positions of 0. 

As a specific illustration of a solution of the above form one might consider what 
happens when any type of bottleneck, a slow car, a road block, traffic light, etc. is 
suddenly removed at time 0. Suppose, for example, that a traffic light has been red for 
a sufficiently long period of time that long queues have developed in both lanes, and that 
a steady flow of traffic is entering the highway from the side roads. If at time zero, the 
light turns green, the initial conditions on the left of the stop-light will be represented 
by some point along the right-hand boundary of the hodograph plane, and the initial 
conditions to the right of the stop-light will map into a point on the upper boundary. 
The solution for this type of initial data is shown in Fig. 5 and consists of a minus wave 
to the singular characteristic, a wave solution along the singular characteristic and then 
a plus wave, each separated by regions of constant density. This solution will be valid, 
however, only until such time when the light turns red again, the queues of traffle become 
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Fig. 5. The hodograph map shown by the broken line in (a) for 8 = 0.1 and the physical plane 
solution (p and q vs y) in (b) represent the flow of traffic immediately after a stop light turns green. 
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empty, or interactions occur between the waves resulting from successive cycles of the 
traffic light. 

4. Other types of initial data. There are many types of initial data other than 
the discontinuous type described in Sec. 3 that yield relatively simple solutions. We 
consider here a few other examples. We will assume, however, that the initial data are 
such that the solutions lie entirely within the hyperbolic region of the hodograph plane. 
We have not been able to determine what should happen if the initial data lie partly 
in the elliptic region or if the solution touches the boundary of the elliptic region at 
some later time. 

One simple class of solutions known as simple waves [6] arises if the image of the 
initial data in the hodograph plane lie entirely along a single hodograph characteristic. 
In this case, one can show that unless shocks develop the solution maps into this same 
hodograph characteristic for all times and that p and q are constant along the physical 
plane characteristics. The velocity w of the physical plane characteristics, Eq. (4), being 
a function of only p and q will also be constant and to obtain w, one chooses the plus or 
minus sign in Eq. (4) accordingly as the wave lies on a plus or minus hodograph character- 
istic. 
To construct the solution from the initial data p(x, 0), q(x, 0), one draws through 
each point x at time 0, the straight line physical plane characteristic with the velocity 
w determined by the value of p(x, 0), g(x, 0) and assigns the same p and q to all points 
on this characteristic. If none of these physical plane characteristics intersect, the 
complete solution is thereby uniquely determined. If some characteristics do intersect, 
a shock develops at the earliest time this happens and the shock path must, in most 
cases, be determined numerically. 

Suppose now we assume that the initial data do not lie on a hodograph characteristic 
but that p(x, 0) and q(z, 0) are constants for x = 0 and for x = 2» with values (p, , qg,) 
and (p2 , g2) respectively, and that for0 < x < x, the map of the initial data into the 
hodograph plane lies on an are between (p, , g,) and (pz , q2) which is never tangent to a 
hodograph characteristic and on which x(p, q) is single valued. Initial data of this type 
and its image in the hodograph plane are illustrated in Figs. 6(a) and (b). 

To find the solution in this case one proceeds in the following way. From Eq. (3) 
we calculate the values of the partial derivatives, ¢, and ¢, , on the initial data curve 
and from these deduce which side of the initial data curve corresponds to positive times. 
If one of the hodograph characteristics drawn through the point (p, , g,) intersects one 
of the characteristics drawn through (p. , g2) and the intersection, (ps3 , q3), lies on the 
side of the initial data curve corresponding to positive time, these three curves will 
enclose a domain of dependence in the hodograph plane as illustrated in Fig. 6(b). Such 
intersections are not always assured since the characteristics may end at the perimeter 
of the ellipse without intersecting each other. If such a domain does exist, however, 
the values of ¢(p, g) 2 0 and x(p, q) can be determined for all points inside this domain 
using the method of characteristics [10], or some alternative scheme. 

This solution can be mapped into the physical plane by drawing the curves of con- 
stant time in the hodograph plane. Along these curves p and q may be considered as 
functions of x and be plotted as such for various times. The limiting values of x and / 
obtained as we approach the boundary of the hodograph domain of dependence from 
its interior define two curves in the (2, ¢) plane which emanate from the points (0, 0) 
and (z, , 0) and have an intersection at the value of x and ¢ corresponding to the point 
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(ps , 93) in Fig. 6(b). These curves are the physical plane characteristics and define the 
physical plane domain of dependence illustrated by the area Oaz, , region 3, in Fig. 
6(c). The solution in the hodograph domain of dependence determines the solution 
everywhere in the physical plane domain of dependence. 

To find the physical plane solutions outside the domain of dependence, we make 
use of the fact described above that if the map of a physical plane solution coincides 
with a hodograph characteristic, the solution is a simple wave. Since the solutions along 
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Fic. 6. An example of continuous initial data, together with their hodograph and physical plane domains 
of dependence are shown in (a), (b) and (c) respectively. The solutions, p and qg, are constants in regions 
1, 4 and 6 of (c) and simple waves in regions 2 and 5. 


Oa and z,a of Fig. 6(c) map into the hodograph characteristics, (p, , q:) — (Ps , gs) and 
(po , G2) — (Ps , Gs), Tespectively, it is possible to continue the physical plane solution 
to points outside the region 3 by means of simple waves. From any pcint on 0a, for 
example, at which p and q have been found, one can draw the appropriate straight 
line characteristic and assign to every point on this characteristic outside region 3 the 
same values of p and gq. This gives the solution in region 2 of Fig. 6(c). Similarly the 
solution in region 5 is obtained from a simple wave originating from 29a. In regions 1, 4 
and 6, p and q will be constant and have the values (p, , ¢:), (ps , Qs) and (Pe , 2), re- 
spectively. 

Since the solution in the hodograph domain of dependence can be found by means 
of the method of characteristics, we can always obtain an approximate solution for this 
type of initial data at least by numerical means. For certain types of initial data, how- 
ever, one can also obtain the solution analytically. If p(x, 0) and g(x, 0) are linear in z 
for 0 < x < 2, there exist solutions of Eq. (1) having the form 


p(x, t) = a,(t) + b,(dz, q(x, t) = a(t) + b(dz. (13) 
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If one substitutes this into Eq. (1), one obtains the four ordinary differential equations 


b, — 2a,b, — B(a,b. + a.b,) = —da,/dt, (14a) 
b, — 2a.b. — B(a,b. + a,b;) = da,/dt, (14b) 
2b,(b,; + Bb.) = db,/dt, 2b.(b. + Bb,) = —db./dt, (14¢, d) 


for a,(t), b, (2), a.(t) and b,(t) whose values at ¢ = 0 are known from the initial conditions. 
It is even possible to find the complete analytic solution of Eq. (14). If we divide Eq. 
(14c) by (14d), the resulting equation is homogeneous of degree zero and can be inte- 
grated to give a relation between b,(¢) and b.(¢). From this one can then find the solutions 
b,(¢) and 6,(¢) of Eqs. (14¢, d). Substitution of this into Eqs. (14a, b) then gives a pair of 
simultaneous linear differential equations which also can be integrated. 

In most cases, the solutions of Eq. (14) can not be expressed in terms of elementary 
functions, but there are two special cases in which the solutions are very simple. If the 
hodograph map of the initial data is parallel with the singular characteristic, b,(0) = — 
b,(0), or parallel with either axis, b,(0) = O or b,(0) = O, then it remains this way 
b, (4) = — b(t), b, (2) = 0, or b.(t) = 0, respectively, for all time. The complete solutions 
for these cases are given in [8]. 

The above solutions apply only within the physical plane domain of dependence, 
region 3 of Fig. 6(c), or equivalently in the hodograph domain of dependence determined 
by the two characteristics through (p, , g,) and (pz , g2). The solutions outside these 
domains are again found by using simple waves. The solutions for two examples with 
b(t) = 0 are illustrated in Figs. 7 and 8. 

In Fig. 7, waves start to develop at x = 0 and x = 2, at time zero. At time f, , p 
and q each show a region of constant density, a wave, a region of linearly varying density, 
another wave, and finally constant density again in that order as one goes from left to 
right. The waves develop at the expense of the linear region until at time ¢, , the latter 




















Fic. 7. The solution corresponding to the data shown for ¢t = 0 is represented by plotting p and q as 
a function of z for several values of the time t) = 0 < t; :-: <t,. 
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Fic. 8. This figure is the same as Fig. 7 except that the initial data lead to the formation of shocks. 


The final configuration (not shown) eventually becomes two diverging shocks. 


region disappears. At time ¢, , we have two waves separated by a region of constant 
density (ps , g,), and for subsequent times the two waves continue to travel away from 
each other and spread over wider ranges of .r. 

The analogous situation but with p(x) increasing with x at ¢ = 0 is shown in Fig. 8. 
At time /, , the region of linearly varying density has vanished and at time ¢,; we have two 
waves separated by a region of constant density (p; , g;). The two waves, however, are 
of the type that develop shocks and one of them is about to form at time ¢, . Eventually 
one obtains a configuration with just two shocks separated by a region of constant 
density. The constant densities are, however, determined by the intersection of the 
two shock polars for (p, , g,) and (p2 , g,) and are not (ps; , q;) which are determined by 
the intersection of the characteristics through (p, , q,) and (p. , q:). This final state is 
achieved by means of waves traveling back and forth between the two shocks, decreasing 
in amplitude and increasing in width with each trip until the densities are nearly con- 
stant. The details of this are rather tedious to follow but for 8 «< 1 these secondary 
waves will be of very small amplitude and are almost completely absorbed each time 
they overtake a shock. 

5. Acknowledgments. The authors are indebted to Professor W. Prager for suggest- 
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it, and to Professor R. Meyer for many very helpful suggestions. 


REFERNCES 
1. M. J. Lighthill and G. B. Whitham, On kinematic waves. II. A theory of traffite flow on long crowded 
roads, Proc. Roy. Soc. A229, 317-345 (1955) 
2. P. I. Richards, Shock waves on the highway, J. Opns. Research Soc. Amer. 4, 42-51 (1956) 
3. S. C. De, Kinematic wave theory of bottlenecks of varying capacity, Proc. Cambridge Phil. Soc. 52, 


564-572 (1956 








204 J. H. BICK AND G. F. NEWELL [Vol]. XVIII, No.2 


4, K. Ostwatitsch, Uber die Charakteristikenverfahren der Hydrodynamik, Z. angew. Math. Mech. 25, 
195-208 (1947) 

5. P. D. Lax, Non-linear hyperbolic equations, Communs. Pure and Appl. Math. 6, 231-258 (1953) 

6. R. Courant and K. O. Friedrichs, Supersonic flow and shock waves, Interscience Publishers, New 
York, 1948, Chaps. II and IIT 

7. E. Kamke, Differentialgleichungen, Lésungsmethoden und Lésungen, Akademische Verlagsgesell- 
schaft, Leipzig, 1944, p. 32 

8. J. H. Bick and G. F. Newell, A continuum model for traffic flow on an undivided highway, Brown 
University report IBM-22 (mimeographed), 1958 

9. G. Birkhoff, Hydrodynamics, Dover Publications, New York, 1955, Chap. IV 

10. I. Petrovsky, Lectures on partial differential equations, Interscience Publishers, New York, 1954, 


p. 68 





—NOTES— 


UPPER BOUNDS AND SAINT-VENANT’S PRINCIPLE IN 
TRANSIENT HEAT CONDUCTION* 


By BRUNO A. BOLEY (Institute of Flight Structures, Columbia University, N. Y.) 


Summary. An investigation is carried out on transient heat-conduction problems 
with prescribed surface temperature, and the validity of Saint-Venant’s principle in 
parabolic boundary-value problems is discussed. 

Some upper bounds for the steady-state temperature in a body whose surface temper- 
ature is prescribed over a small portion of the bounding surface were derived in a recent 
investigation [1] concerned with various possible formulations of Saint-Venant’s principle. 
It was indicated there that this principle was a general property of elliptic differential 
equations, and would not hold all the time for differential equations of other types. For 
the parabolic type, however, it appeared that such a principle could be valid in a rather 
general form; this conclusion was borne out by earlier work [2] in which actual use of 
the principle was made. It is the purpose of this work to investigate this matter a little 
further, and to establish for the transient heat conduction case (which is of the parabolic 
type) the same sort of upper bounds which were treated in [1]. 

We note first a number of properties of fundamental solutions of the Fourier heat- 


conduction equation 


KV*7T(Q, t) = wnt: (1) 


which follow from the fact that [3]** 
if 7T(P,t)>0 and 7(Q,0) = 0, then 7(Q,t) > 0 (2) 


or in other words that it is impossible to lower the temperature anywhere in the interior 
of a body by raising its boundary value. 

Consider now the fundamental solution corresponding to a unit source liberated at 
Q, at a time 7 and with vanishing temperature on the surface (analogous therefore to a 
Green’s function); the temperature at Q, due to this source is denoted by G(Q; , Q2 , 
t — r), and then 

G(P,Q.,t— 7) = 0; GQ, ,Q, ,0) = 0. (3) 

The following properties can now be derived: 

(a) G(Q, , Q.,t — r) > 0, since Eqs. (3) hold, and G is positive by definition in the 
neighborhood of Q, = Q. , and regular everywhere but at Q, = Q, ; therefore Eqs. (2) 
apply and give immediately the desired result. 





*Received February 16, 1959. This work is part of a project supported by the Office of Naval Re- 


search. 
**We will always use P to denote a point on the boundary, and Q a generic point, either on the 


boundary or in the interior. 
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(b) It follows from (a) that 


aG(P, Q. , t — 7) 


> 0 
On, ~ 


where n, is the inward normal to the surface at P. 

(c) Consider two domains, D and D*, where D is completely contained in D*, and 
let G and G* be the corresponding Green’s functions; then G*(Q, , Q. ,4 — r) > GQ. , 
Q. , 4 — r) with Q, and Q, in D. The difference (G* — G) is in fact initially zero for Q, 
in D, and its value on the surface of D is G* because of (a); hence Eqs. (3) apply and the 
desired result follows. 

(d) Consider now the case in which the surfaces of the two domains just described 
have a common portion; it then follows from (c) that on the common boundary 


dG*(P, Q. , t — 7) ‘ dG(P, 7 ,t— tT) 
on, — on, 


Note that all the above properties are well-known in the steady state case [4]. 
The solution due to a prescribed surface temperature is [5]: 


me r i IG(P,Q,t — : 
T(Q,t) =« | | TP, 1) Se” asp) dr, (4) 
Jo ds On, 
where S is the surface of the body, and where it has been assumed that 7(Q, 0) = 0. 


We wish to consider a problem in which 7'(P, t) is zero everywhere except in a (small) 
portion S, , where we may set 


tr, Oi S&S FhH ST : (5) 


Hence with property (b) we have 


aG(P, Q, t — 7) 


On, 


T(Q,t)|<« | T(r) | dS(P) dr (6) 


J0 J So 


and, further, if D and D* are two domains related as in (d) above: 


* > 6G*(P,Q,t— 7) 5. 
ith 2 d J aS(P) dr (7) 


10,0) |S | T’(r) | - 
Jo JSo on, 

Eqs. (6) and (7) represent upper bounds to the solution. The determination of G or G* is, 

however, often very difficult [6], but, as in [1], it may be at times possible to take for D* 

a half-space (say x > 0 in Cartesian coordinates); then 


1 f R? R® 7 
G*(Q, , Q. ,t — 7) = =— sas exp | -— — exp| —; 8 
ws 5 @ i S[ax(t — 7r)]’~ exp | 4k(t — 5 | exp | A(t — 5} 8) 


where F is the distance between Q, and Q. , and R* the distance between Q, ‘and the 
reflection of Q, in the plane x = 0. With the aid of (8) the integrations in Eqs. (6) and 


2 
Z = 


(7) can be carried out; for example, if S, is the surface bounded by the circle y° + z 
2 - 22 . ee. 
h’/x, whose area is h’, then the upper bound along the line y = z = 0 is given by 


1.0 <et{akee]-(i+8)"et[ado (+S) « 
T. < erf (xt) | ~ an. erf Gxt)? Leo (9) 
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where the constant 7’, has been used in place of the function 7,(¢) so as to avoid the 
choice of specific applied-temperature history. The right-hand side of (9), a function of 
the two parameters (2/h) and [h/(4 « t)'], is plotted in Fig. 1. Note that ast > » the steady 
value for the upper bound is reached, i.e. 1 — [1 + (h?/xx*)]~! as found in Ref. [2]*. 
For (x/h) = 0 the upper bound is always unity, and as (1/h) — © it always approaches 
zero; this approach is quite rapid for all values of the parameter h/(4 « ¢)', so that the 
temperature for (x/h) > 1 could possibly be considered small compared to 7, . In this 
case we may then say that Saint-Venant’s Principle holds, and it is of interest to note 
that it holds more strongly for the transient than for the steady state, or, in other words 
that the steady state values provide an upper bound for the transient cases, in agreement 
with [2]. This type of conclusion cannot be extended to problems in which the transient 
case is governed by a hyperbolic equation [7]. 
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